LoginSignup
1
4

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

Posted at

概要

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

自己相関係数

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) と呼びます.

image.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$を作る
  2. $x_{t-1},\cdots,x_{t-j+1}$から$x_{t-j}$を予測する重回帰モデル$f_2$を作る
  3. 両モデル$f_1,f_2$の残差$\boldsymbol{e}_1,\boldsymbol{e}_2$を計算する
  4. 残差の相関係数$\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) といいます.

image.png

ライブラリを使った計算

実は,上記のような実装をしなくても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()

image.png

ソースコード

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

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

1
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
4