概要
シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、自己相関係数・偏自己相関係数について扱います。
自己相関係数
1次元の離散的な時系列データ$\{x_t\}_{t=0}^{\infty}$の $j$次の自己相関係数（auto correlation coefficient） とは，次の2つの系列データの相関係数のことです．
- 元の時系列データ$\{x_t\}_{t=0}^{\infty}$
- 元の時系列データから時間$j$だけずらした（シフトさせた）データ$\{x_{t+1}\}_{t=0}^{\infty}$
例えば，
..., 0, 1, 2, 3, 2, 1, 0, ...
のように「0,1,2」が周期的に並ぶデータの1次の自己相関係数を計算したいなら，
- 元データ「..., 0, 1, 2, 3, 2, 1, 0, ...」
- 元データを1つぶんずらしたデータ「..., 1, 0, 1, 2, 3, 2, 1, ...」
の相関係数を計算すれば良いです．
なお，0次の自己相関係数は必ず1になります．
自己相関係数をPythonで計算する
気象庁の「過去の気象データダウンロード」から取得したデータの自己相関係数（0次～100次）を計算してみましょう．
地点：京都
期間：2022年11月26日～2023年11月26日
データの種類：時別値
前処理方法は割愛します．
変数dfに下表のようなデータが入っていると思ってください．
dfはPandasのDataFrame型です．
| datetime | temp |
|---|---|
| 2022-11-26 01:00:00 | 10.9 |
| 2022-11-26 02:00:00 | 10.9 |
| 2022-11-26 03:00:00 | 10.7 |
| 2022-11-26 04:00:00 | 10.7 |
| 2022-11-26 05:00:00 | 10.6 |
| ... | ... |
時間をずらしたデータを作る
1次～100次の自己相関係数を求めるために，df['temp']を1～100だけずらしたデータdf['lag1']～df['lag100']を作成します．
これを行うには，Pandasの.shift()メソッドが便利です．
MAX_LAG = 100 # 100次まで求める
df_lag = df[['temp']].copy() # dfからカラム「temp」を取り出し，df_lagへ格納
for lag in range(1, MAX_LAG+1):
df_lag[f'lag{lag}'] = df_lag.temp.shift(lag) # カラム「temp」をlagだけ下へずらす
df_lag = df_lag.dropna() # ずらすと先頭にnanが出るのでdropする
これで，1～100個ずらしたデータがまとめてdf_lagへDataFrameとして格納されました．
自己相関係数を求める
あとはdf_lagに対して相関係数を求めてやれば，それが自己相関係数となります．
相関係数は.corr()というPandas DataFrameのメソッドが使えます．
corr = df_lag.corr()
corr = corr.loc[['temp'], :]
求まった自己相関係数を可視化すると下図のようになります．
このように，自己相関係数を次数の関数で表したものを 自己相関関数（ACF; AutoCorrelation Function） といい，そのグラフを コレログラム（Correlogram） と呼びます．
ちょうど24の倍数次で自己相関関数が極大をとっていますね．
これは気温が日周期で変動していることの反映です．
偏自己相関係数
自己相関係数を使うと，時間差$j$をもつデータ$x_t$と$x_{t-j}$の間の相関関係を知ることができますが，場合によっては見せかけの相関が見えてしまう場合があります．
それは，$x_t$と$x_{t-j}$の間にある$x_{t-1},\cdots,x_{t-j+1}$が$x_t$や$x_{t-j}$に影響を与えている，すなわち交絡因子が存在する場合です．
交絡を除去して自己相関を考えたい場合は，交絡因子である$x_{t-1},\cdots,x_{t-j+1}$が$x_t$の影響を取り除いたうえで，$x_t$と$x_{t-j}$の間の相関係数を求める必要があります．
このように，交絡を除去してもとめた自己相関係数のことを 偏自己相関係数（Partial Autocorrelation Coefficient） といいます．
偏自己相関係数では，次の2つの相関係数を計算します．
- $x_{t-1},\cdots,x_{t-j+1}$から"予測できない"$x_t$の情報
- $x_{t-1},\cdots,x_{t-j+1}$から"予測できない"$x_{t-j}$の情報
ここで「予測できない情報」とあいまいに書きましたが，これを具体的に言うと，
重回帰モデル「$x_t=\beta_0 + \beta_1x_{t-1} + \cdots + \beta_{j-1}x_{t-j+1}$」で説明できない部分，つまり残差「$x_t-\hat{x_t}$」
のことです．
偏自己相関係数をPythonで計算する
時間をずらしたデータを作成するところまでは，自己相関係数のときと同じです．
偏自己相関係数を求めるステップとしては，
- $x_{t-1},\cdots,x_{t-j+1}$から$x_t$を予測する重回帰モデル$f_1$を作る
- $x_{t-1},\cdots,x_{t-j+1}$から$x_{t-j}$を予測する重回帰モデル$f_2$を作る
- 両モデル$f_1,f_2$の残差$\boldsymbol{e}_1,\boldsymbol{e}_2$を計算する
- 残差の相関係数$\text{Corr}(\boldsymbol{e}_1,\boldsymbol{e}_2)$を求める
です．
今回も1次～100次までの偏自己相関係数を求めたいので上記ステップを各lagについて行います．
pacf_values = [1., df_lag.corr()['temp']['lag1']]
for lag in range(2, MAX_LAG+1):
# t-1,t-2,...,t-lag+1からtを予測する
lr_forward = LinearRegression()
lr_forward.fit(df_lag[[f'lag{i}' for i in range(1, lag)]], df_lag['temp'])
# t-1,t-2,...,t-lag+1からt-2を予測する
lr_backward = LinearRegression()
lr_backward.fit(df_lag[[f'lag{i}' for i in range(1, lag)]], df_lag[f'lag{lag}'])
# lr_forwardの残差とlr_backwardの残差の相関を求める
residual = pd.DataFrame()
residual['residual1'] = df_lag['temp'] - lr_forward.predict(df_lag[[f'lag{i}' for i in range(1, lag)]])
residual['residual2'] = df_lag[f'lag{lag}'] - lr_backward.predict(df_lag[[f'lag{i}' for i in range(1, lag)]])
# 結果をリストへ格納
pacf_values.append(residual.corr()['residual1']['residual2'])
求まった偏自己相関係数を可視化すると下図のようになります．
このように，偏自己相関係数を次数の関数で表したものを 偏自己相関関数（PACF; Partial AutoCorrelation Function） といいます．
ライブラリを使った計算
実は，上記のような実装をしなくてもstatsmodelsを使えば一発で計算できます．
# ライブラリで自己相関係数，偏自己相関係数を求める（値のみ）
from statsmodels.tsa.stattools import pacf, acf
pacf_values = pacf(df['temp'], nlags=MAX_LAG)
acf_values = acf(df['temp'], nlags=MAX_LAG)
おまけにグラフを描く関数も用意されています．
# ライブラリで自己相関係数，偏自己相関係数を求める（グラフ付き）
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
plot_acf(df['temp'], lags=MAX_LAG, ax=ax1)
plt.xticks(range(0, MAX_LAG+1, 6))
plt.grid()
ax2 = fig.add_subplot(212)
plot_pacf(df['temp'], lags=MAX_LAG, ax=ax2)
plt.xticks(range(0, MAX_LAG+1, 6))
plt.grid()
plt.show()
ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．
本記事のコードは01-JMA_data_analyticsの中に入っています．