【気象データで時系列解析②】自己相関係数・偏自己相関係数

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

# 自己相関係数
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で計算する

気象庁の「[過去の気象データダウンロード](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)」から取得したデータの自己相関係数（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()`メソッドが便利です．

```python
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のメソッドが使えます．

```python
corr = df_lag.corr()
corr = corr.loc[['temp'], :]
```

求まった自己相関係数を可視化すると下図のようになります．
このように，自己相関係数を次数の関数で表したものを **自己相関関数（ACF; AutoCorrelation Function）** といい，そのグラフを **コレログラム（Correlogram）** と呼びます．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/7466e259-4959-0f27-a734-2869f4a822d9.png)

ちょうど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で計算する

時間をずらしたデータを作成するところまでは，[自己相関係数のとき](#時間をずらしたデータを作る)と同じです．

偏自己相関係数を求めるステップとしては，

1. $x_{t-1},\cdots,x_{t-j+1}$から$x_t$を予測する重回帰モデル$f_1$を作る
1. $x_{t-1},\cdots,x_{t-j+1}$から$x_{t-j}$を予測する重回帰モデル$f_2$を作る
1. 両モデル$f_1,f_2$の残差$\boldsymbol{e}_1,\boldsymbol{e}_2$を計算する
1. 残差の相関係数$\text{Corr}(\boldsymbol{e}_1,\boldsymbol{e}_2)$を求める

です．

今回も1次～100次までの偏自己相関係数を求めたいので上記ステップを各lagについて行います．

```python
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）** といいます．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/ea130de5-6e94-64a0-e446-db1fdd0dda27.png)

# ライブラリを使った計算

実は，上記のような実装をしなくても`statsmodels`を使えば一発で計算できます．

```python
# ライブラリで自己相関係数，偏自己相関係数を求める（値のみ）
from statsmodels.tsa.stattools import pacf, acf

pacf_values = pacf(df['temp'], nlags=MAX_LAG)
acf_values = acf(df['temp'], nlags=MAX_LAG)
```

おまけにグラフを描く関数も用意されています．

```python
# ライブラリで自己相関係数，偏自己相関係数を求める（グラフ付き）
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()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/ba3af8d9-3df1-395b-bc5e-1e7c2b256700.png)


# ソースコード

Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．

本記事のコードは`01-JMA_data_analytics`の中に入っています．

https://github.com/matsuda-tkm/qiita-calender-2023