概要

シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、自己相関係数・偏自己相関係数について扱います。

自己相関係数

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の中に入っています．