非定常データにはccf()を使わない方がいい

結論
ccf()は全てのlagで同一の平均・サンプル数を用いて相関係数を算出するため、定常過程であることを前提にしている。
非定常データで相互相関を見る場合はccf()を使わず、lagごとに平均とサンプル数を出してcorr()した方がよい。

設定
x（先行要因）y（結果）についてt=1～30のデータがあるとします。
相関のありそうな適当な非定常時系列データを作ります。

python

np.random.seed(0)
y = np.random.randn(30)
z = np.random.randn(30)
x = 1.8 * y + np.sqrt(1 - 0.8**2) * z
idx = np.argsort(y)
y = y[idx]
x = x[idx]

相互相関出してみる
statsmodelsライブラリのccf()だと
※引数がccf(遅行、先行)の順番になることに注意。

python

ccf(y,x)[0:6]

python

array([0.96105211, 0.75651025, 0.67512154, 0.59980774, 0.49731567, 0.40271129])

一方同じデータでshift()corr()すると

python

df = pd.DataFrame({'x':x, 'y_lag0':y})
df['y_lag1'] = df['y_lag0'].shift(-1)
df['y_lag2'] = df['y_lag0'].shift(-2)
df['y_lag3'] = df['y_lag0'].shift(-3)
df['y_lag4'] = df['y_lag0'].shift(-4)
df['y_lag5'] = df['y_lag0'].shift(-5)
df.corr()['x']

※ccf()ではcorr(y[t+k], x[t])を計算してくれるが、これは「t時点のxがt+k時点のyに及ぼす影響」を評価している。（遅行側を過去方向にシフトしている）。shift()でこれに合わせるにはyを上側にずらすので、y.shift(-k)となる。（shift(k)だと下方向にずらしてしまう）

python

x         1.000000
y_lag0    0.961052
y_lag1    0.912261
y_lag2    0.899586
y_lag3    0.886459
y_lag4    0.849284
y_lag5    0.798564
Name: x, dtype: float64

まとめるとlag1以降全然違うことがわかる

	ccf()	corr()
lag0	0.96105211	0.961052
lag1	0.75651025	0.912261
lag2	0.67512154	0.899586
lag3	0.59980774	0.886459
lag4	0.49731567	0.849284
lag5	0.40271129	0.798564

どうして一致しないのか
ccf()とshift()corr()では計算に使うサンプルサイズと平均・分散が違う。
相関係数の定義は以下であるが、

ccf(k)=\frac{\sum_{t=0}^{N-k-1}(x_{t+k}-\bar x)(y_t-\bar y)}{N\cdot\sigma_x\sigma_y}

ccf()はlag0 の全サンプルから計算されたパラメータを全lagの計算で再利用している。
上記のうちN（サンプル数）、$\bar x$・$\bar y$（平均）、$\sigma_x$・$\sigma_y$（分散）は全サンプルから算出された（lag0時点の）値を全てのlagで使い回している。
上記の例だとccf()ではlagがいくつであろうとN=30固定である。
分母積和の方は行数で変わるため、分母が固定されているとccf自体はlagが進むほど小さな値になる。
一方shift()corr()では（欠損行は勝手に除かれるので）lag1ならN=29、lag2ならN=28を使う。

またccf()では上記定義のとおり平均・分散も固定値を再利用するため、単位根のあるデータを使うとトレンドによっては過小評価or過大評価することになる。
shift()corr()では当然lagごとに都度平均を算出し直すので、確実に相関係数に関わったデータのみのパラメータが使用されている。

使い分け
サンプルサイズが大きく、定常性があるデータならどっちを使ってもほとんど変わりないはず。
ccf()がいい例
・単純に探索目的
・何らかの変換をして定常性を確保できるなら（ただし単純に差分を取るのはトレンド除去と同時に蓄積効果や低周波成分も消し飛ばすことに注意。「変化量の同期」を見たいならいいが）
shift()corr()がいい例
・定常性のないデータを使いたいとき
・サンプルサイズが小さいとき
・そのlagで実際に相関があるか議論したいとき

以下実験（見なくていい）
・3ラグ遅れて影響が出るデータ
実例考えたけど太陽のフレア観測からある数値データ観測までのラグ検出くらいしか思い浮かばなかった。
検出ラグの期間を推定したいみたいな...

python

np.random.seed(0)
n = 1000
t = np.arange(n)
noise = np.random.randn(n)
x = np.random.randn(n)
y = np.roll(x, 3) * 0.8 + 0.3 * np.random.randn(n)

	ccf()	corr()
lag0	0.00313441	0.003134
lag1	-0.02748165	-0.027483
lag2	-0.02262951	-0.022662
lag3	0.93891773	0.939805
lag4	-0.02980801	-0.029841
lag5	-0.01283555	-0.012842

ちゃんとlag3を検出できる。

・「変化量」が影響するデータ
金利と株価とかがいい例なんじゃないかな
投資家の心理として金利が〇%になったから株買う、ではなく金利下がったから買う（変化後の数値ではなく変化量の正負に応答している）みたいな

python

n = 1000
t = np.arange(n)
noise = np.random.randn(n)
dx = np.random.randn(n)
dy = dx * 0.8 + 0.4 * np.random.randn(n)
x = np.cumsum(dx)
y = np.cumsum(dy)

↓普通の相関係数

	ccf()	corr()
lag0	0.83270227	0.832702
lag1	0.82367074	0.826603
lag2	0.81390974	0.819863
lag3	0.80463824	0.813273
lag4	0.79599233	0.807226
lag5	0.78781568	0.801314

↓差分取った相関係数。明らかにdiff1に相関していることがわかる

	ccf()	corr()
lag0	0.89511438	0.895114
lag1	0.05145265	0.051409
lag2	-0.01258011	-0.012570
lag3	-0.05118882	-0.051098
lag4	-0.00364812	-0.003644
lag5	0.03998239	0.039899

・蓄積効果があるデータ
薬投与量と何らかの観測値とか、じわじわ効いてくる系
相互相関で見るもんじゃない

python

n = 1000
t = np.arange(n)
noise = np.random.randn(n)
x = np.random.randn(n)
kernel = np.exp(-np.arange(10) / 3)
kernel /= kernel.sum()
y = np.convolve(x, kernel, mode="same") + 0.2 * np.random.randn(n)

↓普通の相関係数。どう見ても相互相関がありそうには見えない...。

	ccf()	corr()
lag0	0.15081579	0.150816
lag1	0.0832125	0.083178
lag2	0.04179009	0.041848
lag3	0.04429092	0.044359
lag4	0.03500977	0.035046
lag5	0.0483471	0.048444

↓差分取った相関係数。ただでさえ薄く広くな低周波成分が完全に消えた

	ccf()	corr()
lag0	-0.04620924	-0.046209
lag1	-0.03095754	-0.030934
lag2	0.0275225	0.027489
lag3	0.01265399	0.012626
lag4	-0.04838557	-0.048252
lag5	0.02766006	0.027621