【実践時系列解析】電力需要をVARモデルで予測する

概要

VARモデルとはARモデルの多変量版のこと。
数式で表すと

y_{t} = c + \sum_{p=1}^{N} \varphi_{p} * y_{t-p} + \varepsilon_{n}

c：定数項
$\varphi$：q✖️qの回帰係数の行列
$y_{t-p}$：q✖️１次のベクトル行列。qは使用する時系列の数で決まる
$\varepsilon$：誤差項

このモデルを利用すれば、目的の時系列以外の時系列データをモデルに組み込むことができ、予測や因果推論することができる。

当記事はVARモデルを適用のための手順を、時系列の可視化からモデルの評価まで行う

データセット

Prediction of electric power
Let's share code and approaches so as to improve our data-science skill each other.

import

from statsmodels.tsa.api import VAR
from statsmodels.stats.stattools import durbin_watson
from statsmodels.tsa.stattools import adfuller
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 階差を戻す作業、VAR
def invert_forecast_data(forecast_data, df_train):
    tmp_df = forecast_data.copy()
    col = df_train.columns
    for i in col:
        tmp_df[i] = df_train[i][-1] + forecast_data[i].cumsum()
    return tmp_df

前処理

train = pd.read_csv('data/train_df.csv')
train['DATE'] = pd.to_datetime(train.DATE)
group_train = train.groupby('DATE').mean()
group_train.sort_values(by='DATE', inplace=True)
group_train

使用変数の選定

とりあえず全ての時系列を可視化する

col = group_train.columns
fig, axes = plt.subplots(len(col), 1, figsize=(10, 4*len(col)), sharex=True)

for i, column in enumerate(col):
    sns.lineplot(data=group_train, x='DATE', y=column, ax=axes[i])
    axes[i].set_ylabel(column)

plt.xlabel('DATE')
plt.show()

必要そうな変数を調べるために各時系列の相関を確認する

sns.heatmap(group_train.corr())

目的変数である”POWER”と相関が高いのはV＿2とV＿5。

階差をとったものとも相関を確認。

# 階差をとったものとの相関を確認
col = group_train.columns
for i in col:
    print(i)
    print(group_train[i].shift(-1).corr(group_train.POWER))

#ID
#-0.07014896378483501
#V_1
#0.12822735771878305
#V_2
#0.5921549145887411
#V_3
#-0.4611709346475451
#・
#・
#・

今回は演習なので、とりあえずV＿2とV＿5、POWERのみを使ったVARモデルを構築する

ADF検定で定常性を確認

ARモデルを構築する際の過程に、その時系列が定常過程であるというものがある。故に、ARモデルの多変量版であるVARでも定常性を確認する必要がある。定常性の確認の検定はADF検定。

col = ['V_2','V_5','POWER']
# ADF検定、時系列に何の前処理も施していない場合
for i in col:
    print(adfuller(group_train[i])[1])
#0.6665057151336248
#0.46780545096786963
#0.00016847189974906896

何も前処理を施さない場合、V＿2とV＿5は帰無仮説（系列データは単位根過程である）を棄却できず、定常過程とは言えない。

単位根過程については下記事参照。一言だとランダムウォークのこと。

１回階差をとったもののADF検定の値を確認。

for i in col:
    print(adfuller(group_train[i].diff().dropna())[1])
#5.1788918051229844e-05
#1.0945497819910484e-22
#1.6209693783433674e-22

全てp値は低く、帰無仮説を棄却でき、時系列は定常過程であると言える。どうやら対数まで取る必要はなさそうだ。

階差を考慮したモデルを構築

使用するデータセットを訓練とテストに分割

# 訓練データとテストデータを準備
len_test = 24
df_train = group_train.iloc[:len(group_train) -len_test,:]
df_test = group_train.iloc[len(group_train)-len_test:,:]
# 使用変数のみに絞る
df_train = df_train[col]
df_test = df_test[col]

モデルを構築。AIC（情報量の一種）を基準に、これが最小になる次数を自動で見つけ、モデルを構築する

# モデルのインスタンス生成
var_model = VAR(df_train.diff().dropna())
results = var_model.fit(maxlags=40, ic='aic')
print(results.k_ar)
#6

予測

# 差分の時系列を予測
forecast = results.forecast(y=results.endog,steps=len_test)
forecast_data = pd.DataFrame(forecast,index=df_test.index,columns=col)
# 階差を戻す作業
forecast_data = invert_forecast_data(forecast_data,df_train)

fig, axes = plt.subplots(len(col), 1, figsize=(10, 4*len(col)), sharex=True)

for num, i in enumerate(col):
    sns.lineplot(data=group_train, x='DATE', y=i, ax=axes[num])
    sns.lineplot(data=forecast_data, x = forecast_data.index, y = i, c='r',ax=axes[num], label='predicted')
    axes[num].set_ylabel(i)

plt.xlabel('DATE')
plt.show()

あまり予測がうまくいっていない。

次数を４０にしてモデルを構築する。当然過学習が起こり、よりモデルにフィットした予測を出力するはずである

# モデルのインスタンス生成
var_model = VAR(df_train.diff().dropna())
results = var_model.fit(40)

POWERは特にうまく予測できていることがわかる

モデルの評価（偏自己相関の可視化、ダービン・ワトソン検定）

モデルの残差の偏自己相関を確認してモデルが問題ないか確認する。

# 残差の偏自己相関を確認
for i in col:
    sm.graphics.tsa.plot_pacf(results.resid[i], lags=50)
    plt.xlabel('lags')
    plt.ylabel('pacf')

実際に数値でも確認する。DW検定を利用する。

# ダービン・ワトソン統計量DWは、２に近い時自己相関がないと判断されます。２より十分小さい時は正の自己相関、２より十分大きい時は負の自己相関があると判断されます。
for i in col:
    print(durbin_watson(results.resid[i]))
#2.022688712972664
#2.012947749200144
#1.962888388324888

参考文献