Pythonを用いた一般ガウス状態空間モデルによる時系列分析【外因性と季節性を考慮した実装例】

この記事でやること

• Pythonを用いて一般ガウス状態空間モデルの時系列分析の実装方法を記します
• statsmodelsを用いて実装します
• 季節性，外因性を考慮したときの実装を行います
• 理論的な部分は参考記事を紹介するにとどめます
• 詳細な検討は行いません(行えない...)
• 「時系列分析と状態空間モデルの基礎 RとStanで学ぶ理論と実装」
  読んだけどRわかんなくて困った！pythonで書いてくれ！！っていう人

はじめに

時系列分析を**「時系列分析と状態空間モデルの基礎 RとStanで学ぶ理論と実装」** (いわゆる隼本)を用いて勉強中です。
しかし，自分はPythonに慣れているのでRでの実装は飛ばして，本の著者のブログの下記記事を参考にして実装していました。

• Pythonによる状態空間モデル:https://logics-of-blue.com/python-state-space-models/

この記事でPythonによる一般ガウス状態空間モデルの実装方法は把握できたのですが，外因性があるときの実装方法と効果がのっていませんでした。
そこで，別データを使用して外因性があるデータにおける実装と効果についてまとめます。

*状態空間モデルについては同じく本の著者の下記記事をご参照ください。
状態空間モデルの推定方法の分類:https://logics-of-blue.com/python-state-space-models/
カルマンフィルタの考え方:https://logics-of-blue.com/kalman-filter-concept/

使用したライブラリはこちらから。実装はgoogle colabで行っています。

import pandas as pd
import numpy as np
import os
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.pylab import rcParams
import statsmodels.api as sm

##データ準備

• 東京電力の2019年電力需要データを使用:http://www.tepco.co.jp/forecast/index-j.html

df = pd.read_csv("data/juyo-2019.csv",encoding="shift-jis")
df.head()

データが一時間ごとになっているので，一日ごとのデータにまとめ直します。
groupby関数をつかえばすぐにできます。ついでにデータのインデックスを日付に修正します。

df["date"] = pd.to_datetime(df["date"],format="%Y/%m/%d")
df["power"] = df["power"].values.astype(int)
daily_data = df.groupby("date")["power"].sum().reset_index()
daily_data = daily_data.set_index(["date"])
rcParams["figure.figsize"] = [12,5]
daily_data.plot()

データをみると，周期性のあるギザギザがあることがわかります。
これは7日ごとに変動しているので値曜日に関連する周期であろうと推測できます。
つまり休みの日かどうかで電力需要が変化するということですね。

また，曜日だけでなく，祝日かどうかも電力需要に影響を与えることが推測されます。
祝日は周期的に発生していないので，外因性としてモデル化できます。

**まとめると，下記をモデル化することによってより精度よくモデル化しよう，ということになります。

• 周期性:曜日で7日ごとに発生
• 外因性:祝日でランダムに発生**

##季節性，外因性を考慮しない場合

まずは，季節性，すなわち曜日を考慮しないでモデル化してみます。
季節性，外因性を考慮しない，ローカル線形モデルを用いてモデル化します。このモデル化の場合，

#ローカル線形トレンドモデル
mod_trend = sm.tsa.UnobservedComponents(daily_data,"local linear trend")
res_trend = mod_trend.fit(method="bfgs")

rcParams["figure.figsize"] = 15,15
fig = res_trend.plot_components()

図は一番上が実データ，二番目が状態，三番目がトレンドデータを表しています。
周期性を考慮しないとトレンドがぐちゃぐちゃに変形してしまっていることがわかります。

##季節性を考慮した場合

次に周期性を考慮した場合。UnobservedComponentsにseasonalの引数をいれるだけで完了です。
ここでは曜日を考慮するので引数は7としています。
図の四番目に季節性を表すグラフが追加されています。周期的に変動しているのがわかりますね。

#季節性を入れた，ローカル線形トレンドモデル
mod_season_trend = sm.tsa.UnobservedComponents(daily_data.power,
                                              "local linear trend",
                                              seasonal = 7)

res_season_trend = mod_season_trend.fit(
    method='bfgs', 
    maxiter=1000, 
)

rcParams["figure.figsize"] = 15,20
fig = res_season_trend.plot_components()

先ほどとの大きな違いはtrendデータのギザギザがなくなっている点です。
周期性が取り除かれることによってトレンドが先程より精度よくがみえています。
(周期性のグラフの1月あたりの最初のデータが高い値になってしまっているのが気になりますが...)

また，5月のはじめにlevelのグラフが大きく減っていることがわかります。
これは状態が変わっているというよりは，GWの影響だと考えられます。
そこで次は祝日を外因性として組み込んだ場合のモデルを考えてみます。

##季節性，外因性を考慮した場合

最後に外因性も考慮した場合のモデルです。
外因性として祝日を考慮します。祝日フラグをいれるためにjpholidayをimportします。
(ここではgoogle colabを用いた場合のコードです)
に祝日の場合1となるholiday列を追加します。

土曜日，日曜日は元から休みなので，祝日フラグはたてません。

#祝日フラグを入れる，google colabの場合
!pip install jpholiday
import jpholiday

#祝日フラグ
daily_data["holiday"] = 0
for i in range(len(daily_data)):
  if daily_data.index[i].dayofweek != 6: #日曜日でない場合のみ
    date_i = daily_data.index[i] #日付
    if jpholiday.is_holiday(date_i): #祝日の場合
      daily_data.loc[date_i,"holiday"] = 1
daily_data.head()

祝日だけだと連休が考慮できていないですね。
12～1月の冬休み，5月のGW，8月のお盆休みを手動で1に書き換えます。

list1 = daily_data.query("20190101 <= index <= 20190104").index
list2 = daily_data.query("20190429 <= index <= 20190503").index
list3 = daily_data.query("20190812 <= index <= 20190816").index
list4 = daily_data.query("20191230 <= index <= 20191231").index

for i in list1:
  daily_data.loc[i,"holiday"]=1
for i in list2:
  daily_data.loc[i,"holiday"]=1
for i in list3:
  daily_data.loc[i,"holiday"]=1
for i in list4:
  daily_data.loc[i,"holiday"]=1

最後に，外因性をいれたローカル線形トレンドモデルを作成し，プロットしてみます。
引数にexogを入れて，先程作成した祝日フラグを指定します。

#季節性，外因性を入れたローカル線形トレンドモデル
mod_season_exog_trend = sm.tsa.UnobservedComponents(
    daily_data.power,
    "local linear trend",
    seasonal=7,
    exog=daily_data.holiday
)

res_season_exog_trend = mod_season_exog_trend.fit(
    method="bfgs",
    maxiter=1000, 
)

# 推定された状態・トレンド・季節の影響の描画
rcParams['figure.figsize'] = 15, 20
fig = res_season_exog_trend.plot_components()

5月初旬に着目すると，周期性のみの場合に比べて低下がすこーしだけ抑えられていますね。
しかし，祝日フラグを入れたのみではこの低下は完全には抑えられていないようです。

さらなる改善が必要なようですが，今の自分の実力ではこれが限界でした。。。

##まとめ

• statsmodelを用いて季節性，外因性のあるモデルを実装
• 電力需要を対象に，季節性を曜日，外因性を祝日として実装した
• 周期性のモデルは効果が高く実装できた
• 外因性のモデルは祝日だけでは不十分でさらなる改善が必要

最後までお読み頂きありがとうございました。

参考記事

• Pythonによる状態空間モデル:https://logics-of-blue.com/python-state-space-models/
• さくっとトレンド抽出: Pythonのstatsmodelsで時系列分析入門:https://data.gunosy.io/entry/statsmodel_trend