相互相関

05.1.4

相互相関

最終更新 2020-03-11 読了時間 5 分

まとめ

相互相関（cross-correlation）は、2 本の時系列の片方を 1 期ずつずらしながら相関を計算し、「どれだけずらしたときに最も似るか」を調べる手法です。このページでは FRED から鉄・鋼の生産者物価指数 WPU101 と鉄鉱石の生産者物価指数 PCU2122121221 を取り、海運会社 Genco Shipping（GNK）の月次株価と numpy.correlate(mode="full") で相互相関を計算します。相関が最大になる位置 corr.argmax() から系列長を引いた値がラグで、0 なら同時に動き、正なら株価が物価指数に遅れて動くと読めます。計算の前に ax.twinx() で 2 本の線を単位の違う 2 軸に重ね、目で見た印象と数値のラグを突き合わせます。

直感 #

鉄鉱石や鋼材の値段が上がると、それを運ぶばら積み船の需要が増え、海運会社 GNK の株価も上がりそうです。問題は「同じ月か、数か月遅れか」です。相互相関は片方の系列を 1 か月、2 か月とずらして重ね、重なりの良さ（積和）を全部のずらし幅について計算します。積和が最大になったずらし幅がラグで、ラグが 3 なら「物価指数が動いてから 3 か月後に株価が動く」と解釈できます。平均を引いてから計算するのは、両方の系列が常に正の値だと「ずらさなくても全部正×正で積和が大きい」という見かけの相関が出てしまうためです。

2 本の系列を同じ月次・同じ期間に揃えるそれぞれ平均を引く np.correlate(a, b, “full”) argmax の位置 − (len(b) − 1) = ラグラグ 0 なら同時、正なら a が b に遅れる

FRED からデータを取り出す準備 #

生産者物価指数は FRED から取得します。関数 get_fred_data() は FRED データベースのページと同じもので、系列 ID を渡すと日付索引の DataFrame を返します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import japanize_matplotlib
import numpy as np
import os
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from full_fred.fred import Fred

# FRED_API_KEY = os.getenv('FRED_API_KEY')
fred = Fred()
print(f"FRED APIキーが環境変数に設定されている：{fred.env_api_key_found()}")

def get_fred_data(name, start="2013-01-01", end=""):
    df = fred.get_series_df(name)[["date", "value"]].copy()
    df["date"] = pd.to_datetime(df["date"])
    df["value"] = pd.to_numeric(df["value"], errors="coerce")
    df = df.set_index("date")

    if end == "":
        df = df.loc[f"{start}":]
    else:
        df = df.loc[f"{start}":f"{end}"]

    return df

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目日本語ラベルのために japanize_matplotlib を読む
後の図では凡例に「鉄・鋼の生産者物価指数」と日本語を使います。japanize_matplotlib を読み込んで japanize() を呼ばないと、matplotlib の既定フォントに日本語が無く凡例が豆腐（□）になります。
10-11 行目API キーが環境変数にあるかを先に確かめる
Fred() は環境変数 FRED_API_KEY を自動で読みます。env_api_key_found() が False なら以降の取得が認証エラーになるので、最初に True を見ておくと原因の切り分けが早くなります。
15-16 行目文字列で届く日付と値を数値にする
FRED の API は値を文字列で返し、欠測を "." で表します。pd.to_datetime と pd.to_numeric(errors="coerce") で日付型・数値型に変換し、"." は NaN にします。この変換を省くと np.correlate が文字列を受け取って TypeError を出します。
17-22 行目日付索引にして期間を文字列で切り出す
set_index("date") で日付を索引にすると、df.loc["2018-01-01":"2022-05-01"] のように期間を切り出せます。相互相関では 2 本の系列の長さと期間を揃えることが前提なので、この start / end が後で重要になります。

FRED APIキーが環境変数に設定されている：True

鉄・鋼の生産者物価指数（WPU101） #

ソース： https://fred.stlouisfed.org/series/WPU101

Producer Price Index by Commodity: Metals and Metal Products: Iron and Steel は鉄及び鋼の生産者物価指数（月次）です。2018 年 1 月から 2022 年 5 月までの 53 か月分を取得します。

1
2
df_WPU101 = get_fred_data("WPU101", start="2018-01-01", end="2022-05-01")
df_WPU101.head(5)

	value
date
2018-01-01	217.6
2018-02-01	222.1
2018-03-01	227.5
2018-04-01	234.1
2018-05-01	239.0

1
2
sns.set(rc={"figure.figsize": (15, 8)})
sns.lineplot(data=df_WPU101, x="date", y="value")

<AxesSubplot:xlabel='date', ylabel='value'>

Producer Price Index by Commodity: Metals and Metal Products… figure

GNK の株価を月次に揃える #

Genco Shipping & Trading Limited (＄GNK) はばら積み船を運航する海運会社です。Yahoo Finance から CSV で保存した日次株価を読み、物価指数と比べられるように月次平均にして月初日に揃えます。

1
2
3
4
5
6
df_GNK = pd.read_csv("./GNK.csv")
df_GNK["Date"] = pd.to_datetime(df_GNK["Date"])
df_GNK = df_GNK.set_index("Date").resample("M").mean()
df_GNK.index = df_GNK.index + pd.DateOffset(1)  # 日付を月初に調整する
df_GNK = df_GNK.loc["2018-01-01":"2022-05-01"]
df_GNK.head()

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-2 行目CSV の日付列を日付型に戻す
read_csv で読んだ Date 列は文字列なので pd.to_datetime で変換します。この変換が無いと次の resample("M") が「索引が日付型でない」というエラーで止まります。
3 行目日次を月次平均に落として頻度を揃える
物価指数は月に 1 点、株価は月に 20 点前後あります。resample("M").mean() で各月の平均を 1 点にまとめ、両系列を同じ 53 点に揃えます。mean() の代わりに last()（月末終値）を使う選択もあり、どちらを使うかでラグが変わることもあります。
4 行目月末日を月初日にずらして索引を一致させる
resample("M") の索引は月末日（2018-01-31 など）になりますが、FRED の索引は月初日です。pd.DateOffset(1) で 1 日足して翌月 1 日に動かし、索引の見た目を揃えます。1 か月ぶんずれて見えますが、np.correlate は索引を見ず位置だけで計算するので結果には影響しません。
5 行目物価指数と同じ期間に切る
.loc["2018-01-01":"2022-05-01"] で WPU101 と同じ範囲にします。長さが違うと np.correlate の結果の中心位置がずれ、delay = argmax - (len - 1) の計算が成り立ちません。

	Open	High	Low	Close	Adj Close	Volume
Date
2018-01-01	12.5575	13.4355	12.1375	12.7650	11.157512	603175.0
2018-02-01	13.9720	14.5420	13.1310	13.7600	12.027212	510600.0
2018-03-01	12.9925	13.8525	12.6520	13.4300	11.738769	357050.0
2018-04-01	15.0100	15.6925	14.4025	14.7575	12.899099	796075.0
2018-05-01	15.5280	16.3530	14.9620	15.9580	13.948420	631500.0

物価指数（200〜400）と株価（10〜25 ドル）は桁が違うので、左右に別の縦軸を持つ 2 軸グラフで重ねます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
japanize_matplotlib.japanize()

fig = plt.figure()
ax1 = fig.add_subplot(111)
sns.lineplot(data=df_WPU101, x="date", y="value", label="鉄・鋼の生産者物価指数")
plt.legend(loc="lower right", fontsize="20")

ax2 = ax1.twinx()
sns.lineplot(data=df_GNK, x="Date", y="Close", label="GNK", color="r")
plt.legend(loc="upper left", fontsize="20")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1 行目日本語の凡例を出すための 1 行
japanize() を呼ぶと matplotlib のフォントが日本語対応のものに切り替わります。省くと「鉄・鋼の生産者物価指数」が □ の列になります。
4-5 行目1 本目は左軸（ax1）に描く
fig.add_subplot(111) で作った ax1 に WPU101 を描きます。seaborn.lineplot は「現在の Axes」に描くので、この時点では左軸の 200〜400 のスケールが使われます。
8-9 行目twinx() で右軸を作り株価を重ねる
ax1.twinx() は横軸を共有し縦軸だけ独立した ax2 を返します。ここに GNK の Close を赤で描くと、右軸に 10〜25 ドルのスケールが付き、桁の違う 2 本を同じ横軸で比べられます。FRED のページで使った「10,000 倍する」方法より正確です。
6, 10 行目凡例を左右に分けて重なりを避ける
plt.legend を各軸の描画直後に呼び、loc="lower right" と "upper left" で置き場所を分けています。1 回だけ呼ぶと直前の軸の凡例しか出ないため 2 回呼んでいます。fontsize="20" は 15×8 インチの図で読める大きさにする指定です。

GNKの価格推移 figure

実行結果の読み方 #

青（物価指数）と赤（GNK）はどちらも 2020 年前半に底を打ち、2021 年に大きく上昇しています。目で見ると赤（GNK）がやや先に上がり始めたようにも、同時にも見えます。この「どちらが先か」を数値で決めるのが次の相互相関です。

相互相関を求める #

相互相関関数は、ふたつの信号、配列（ベクトル）の類似性を確認するために使われる。関数の配列の結果がすべて1であれば相関があり、すべてゼロであれば無相関であり、すべて −1 であれば負の相関がある。しばしば、相関と略されることがあり、相関係数と似ているために混同することがある。 (出典：Wikipedia)

numpy.correlate(a, v, "full") は、v を左端から右端まで 1 つずつずらしながら a との積和を計算し、長さ len(a) + len(v) - 1 の配列を返します。中央（添字 len(v) - 1）がずらし 0 に対応するので、最大値の位置から中央を引けばラグになります。

ドキュメント： numpy.correlate

1
2
3
4
5
6
7
8
df_GNK["Close_norm"] = df_GNK["Close"] - df_GNK["Close"].mean()
df_WPU101["value_norm"] = df_WPU101["value"] - df_WPU101["value"].mean()

corr = np.correlate(df_GNK["Close_norm"], df_WPU101["value_norm"], "full")
delay = corr.argmax() - (len(df_WPU101["value_norm"]) - 1)
print("ラグ： " + str(delay))

plt.plot(corr)

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-2 行目平均を引いて見かけの相関を消す
株価も物価指数も常に正の値なので、そのまま積和を取るとどのずらし幅でも大きな正の値になり、山がぼやけます。各系列から自身の平均を引いて正負が混ざる形にすると、「同じ方向に動く月」だけが積和を押し上げ、ラグの山がはっきり出ます。
4 行目"full" モードで全部のずらし幅を計算する
"full" は 2 本が 1 点だけ重なる状態から完全に通り過ぎるまでを計算するので、結果は 53 + 53 − 1 = 105 要素になります。既定の "valid" だと同じ長さの系列ではずらし 0 の 1 点しか返らず、ラグを探せません。
5 行目argmax から中央を引くとラグになる
corr の添字 52（= 53 − 1）がずらし 0 です。argmax() がそこより右なら正のラグで「GNK が WPU101 に遅れて動く」、左なら負で「GNK が先行する」と読みます。ここでは 0 が返り、同じ月に動いていると分かります。
8 行目相互相関の形を描いて山の鋭さを見る
plt.plot(corr) で 105 点の相互相関を描きます。横軸 52 付近に 1 つの鋭い山があれば「ラグ 0 で強く連動」、山がなだらかなら「ずらしても似たような相関で、ラグの判定は弱い」と解釈します。

ラグ： 0

[<matplotlib.lines.Line2D at 0x13e252170>]

[] figure

実行結果の読み方 #

ラグは 0 で、相互相関の図は中央（横軸 52）を頂点とする 1 つの山になっています。鉄・鋼の物価指数と GNK の株価は月次で見ると同じ月に動いており、「物価指数を見てから株を買う」という時間差は月単位では存在しないことになります。

鉄鉱石の生産者物価指数（PCU2122121221）で同じことを試す #

ソース：https://fred.stlouisfed.org/series/PCU2122121221

Producer Price Index by Industry: Iron Ore Mining は鉄鉱石の生産者物価指数です。鋼材より上流の原料なので、GNK の株価に対して先行するかを同じ手順で確かめます。

1
2
df_PCU2122121221 = get_fred_data("PCU2122121221", start="2018-01-01", end="2022-05-01")
df_PCU2122121221.head(5)

	value
date
2018-01-01	121.0
2018-02-01	122.5
2018-03-01	122.5
2018-04-01	133.8
2018-05-01	138.7

1
2
3
4
5
6
7
8
fig = plt.figure()
ax1 = fig.add_subplot(111)
sns.lineplot(data=df_PCU2122121221, x="date", y="value", label="鉄鉱石の生産者物価指数")
plt.legend(loc="lower right", fontsize="20")

ax2 = ax1.twinx()
sns.lineplot(data=df_GNK, x="Date", y="Close", label="GNK", color="r")
plt.legend(loc="upper left", fontsize="20")

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

3 行目左軸に鉄鉱石指数を描く
データを df_PCU2122121221 に替えた以外は前の 2 軸グラフと同じ構成です。鉄鉱石指数は 2021 年半ばに 300 近くまで上がったあと急落するので、鉄・鋼より山が尖っています。
6-7 行目GNK は右軸に同じ色で重ねる
twinx() で右軸を作り、color="r" で GNK を前の図と同じ赤にしています。色を揃えておくと 2 枚の図を見比べたときに「赤は常に株価」と迷わず読めます。
4, 8 行目凡例の位置は前の図と同じ
左下に物価指数、左上に GNK という配置を踏襲しています。図を並べて比較する前提なら、凡例の位置も固定しておくほうが読みやすくなります。

Producer Price Index by Industry: Iron Ore Mining figure

1
2
3
4
5
6
7
8
9
df_PCU2122121221["value_norm"] = (
    df_PCU2122121221["value"] - df_PCU2122121221["value"].mean()
)

corr = np.correlate(df_GNK["Close_norm"], df_PCU2122121221["value_norm"], "full")
delay = corr.argmax() - (len(df_PCU2122121221["value_norm"]) - 1)
print("ラグ： " + str(delay))

plt.plot(corr)

コードの解説番号をクリックすると、コード内の該当行がハイライトされます。

1-3 行目鉄鉱石指数も平均 0 に揃える
GNK 側の Close_norm は前の節で作ったものを再利用し、鉄鉱石指数だけ新たに平均を引きます。片方だけ平均を引き忘れると、相互相関の山が片側に傾いて誤ったラグが出ます。
5-6 行目同じ式でラグを計算する
np.correlate の第 1 引数を GNK、第 2 引数を物価指数にする順番を前の節と揃えています。順番を入れ替えるとラグの符号が反転するので、「正なら GNK が遅れる」という読み方を保つために固定しています。
9 行目山の形を鉄・鋼のときと比べる
結果はこちらもラグ 0 です。plt.plot(corr) の山の高さは系列の分散に依存するので 2 枚の図の高さを直接比べることはできませんが、山の位置が中央にあることは共通で、原料の鉄鉱石も月次では株価と同時に動いていると分かります。

ラグ： 0

[<matplotlib.lines.Line2D at 0x13e373fa0>]

[] figure

使いどころと落とし穴 #

np.correlate は正規化しないので、返ってくる値は相関係数（−1〜1）ではなく積和です。山の位置（ラグ）は読めますが高さを別のペアと比べることはできず、相関係数の形で欲しいなら両系列を標準偏差でも割るか、pandas の Series.corr を shift(k) しながら繰り返す方法に替えます。株価のように右肩上がりのトレンドを持つ系列は、トレンドだけで相互相関の山が広がります。diff() や pct_change() で変化率にしてから計算するほうが「動きの連動」を正しく捉えられます。月次に落とすと 53 点しかないので、ラグ 1 の差は 1 点のずれで決まることになり、日次や週次で点数を増やしてから確認するのが安全です。

まとめ #

相互相関は 2 本の系列を 1 期ずつずらしながら積和を取り、最大になるずらし幅をラグとして読む手法です。実装は np.correlate(a - a.mean(), b - b.mean(), "full") と argmax() - (len(b) - 1) の 2 行で、前提として resample("M") と .loc[start:end] で頻度と期間を揃えておきます。鉄・鋼と鉄鉱石のどちらの物価指数も GNK の株価とラグ 0 で連動しており、月次では先行関係は見つかりませんでした。トレンドを除いた変化率や、より細かい頻度で試すのが次の一手です。