はじめに

この記事では、Pythonで状態空間時系列モデリングを行う方法を、Rによる実行例の「翻訳」を通して解説します。題材として季節変動のある時系列データを用います。

状態空間時系列モデルでは、季節性やトレンドといった時系列背後の構造を反映させて直感的にわかりやすいモデリングを行うことができます。

しかし、Pythonで状態空間時系列モデリングを行おうことは、簡単ではありません。 Pythonで統計モデリング全般を実行する場合、ライブラリの選択肢はほぼstatsmodels¹一択となります。しかし、このライブラリはインターフェイスが複雑で使いづらいうえにドキュメントも読みにくいのが現状です。また、書籍やWeb情報も数えるほどしか見つかりません。

一方で、Rであれば学習環境がPythonに比べて恵まれています。RにはKFAS², dlm³といった状態空間時系列モデルに特化した定番ライブラリが存在します。また、解説書籍・Web情報もPythonと比べて多く出ています。

そこで今回は、Pythonによる状態空間時系列モデリングのハードルを下げるために、まずはRでモデリングを行い「お手本」を作ります。そしてRによる出力結果を再現しながらPythonでモデリングの学習を行います。具体的には季節変動をもつデータについてRのKFASで予測と成分分解を行い、同一の結果をPythonでstatsmodelsを用いて再現する方法を解説します。記事の構成は以下のようになっています：

状態空間モデルで季節変動のある時系列を表現する方法を解説する
KFASを用いてRで予測・状態推定を行う
- 長期予測
- 季節成分とレベル成分の分離
KFASによる出力結果をstatsmodelsで再現する

なお、この記事は状態空間時系列モデリングについて解説する全5本の記事の3本目となっています。

f:id:tatamiyatamiyatatatamiya:20210212134022p:plain:w500

状態空間モデルの概要や線型ガウスモデルのカルマンフィルタによる扱いについては、前回までの記事を参照してください。

はじめに
- 忙しい人のために
状態空間モデルによる季節変動のある時系列データの扱い方
KFASによるモデル構築
- 基本的なモデル構築方法
  - モデルの定義
  - ハイパーパラメータ推定
- 予測・推定結果
  - 長期予測
  - 成分ごとの推定
statsmodelsによるモデル構築
まとめ
参考資料

忙しい人のために

レベル成分を $\mu_t$ 、季節成分を $\gamma_t$ とし、以下のようなローカルレベル+季節成分モデルを考える：

${ \begin{aligned} y_t &= \mu_t + \gamma_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, \sigma_\varepsilon^2)\\ \mu_{t+1} &= \mu_t + \eta^{(\mu)}_t, \quad \eta^{(\mu)}_t \sim N(0, \sigma_{\mu}^2)\\ \gamma_{t+1} &= - \gamma_t - \gamma_{t-1} - \gamma_{t-2} + \eta^{(\gamma)}_t,\quad \eta^{(\gamma)}_t \sim N(0, \sigma_\gamma^2) \end{aligned} }$

RのライブラリKFASを用いてモデル構築を行えば、以下の予測・推定結果をpredict関数を用いて取得できる：
- 観測値の長期予測: $P(y _ {T+i}|Y _ {1:T})$
- 成分別の推定・予測
  - レベル成分: $P(\mu_t|Y _ {1:T})$
  - 季節成分: $P(\gamma_t|Y _ {1:T})$
  - レベル成分+季節成分: $P(\mu_t + \gamma_t|Y _ {1:T})$
PythonのライブラリstatsmodelsのUnobservedComponentsクラスを用いれば、KFASと同様のモデリング結果を再現できる。
使い勝手の面では、statsmodelsはKFASに劣る
- KFASと異なり結果取得用のインターフェイスは統一されていない
- 一部の予測・推定量の取得方法が複雑

状態空間モデルによる季節変動のある時系列データの扱い方

この章では、線型ガウス状態空間モデルの枠組みで季節変動のある時系列データを扱う方法について紹介します。

まず前回までの記事で扱った線型ガウスの状態空間時系列モデルについて復習を行います。

その次に線型ガウスのモデルとして最もシンプルなローカルレベルモデルを紹介します。

そのうえで、ローカルレベルモデルに季節成分を付加したモデルの紹介を行います。

線型ガウス状態空間モデル（前回の復習）

状態空間モデルでは、観測値の背後にマルコフ遷移する状態を仮定し、データを元に状態の推定・予測を行います。

tatamiya-practice.hatenablog.com

中でも特に重要なのが、線型ガウスの状態空間モデルです。

tatamiya-practice.hatenablog.com

線型ガウスの状態空間時系列モデルは、以下のように表されます。

${ \begin{aligned} y_t &= Z\alpha_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, H)\\ \alpha_{t+1} &= T\alpha_t + R\eta_t, \quad \eta_t \sim N(0, Q) \end{aligned}\tag{1}\label{eq:dlm} }$

$y_t$ は観測値、 $\alpha_t$ は状態を表すベクトルです。 1番目の式は観測方程式、2番目の指揮は状態方程式といいます。

$Z, T, R$ はそれぞれ行列で、モデルの構造により決まるパラメータです。

$\varepsilon_t, \eta_t$ はそれぞれ観測ノイズ、状態ノイズと呼ばれます。両ノイズともに正規分布に従い、共分散行列はそれぞれ $H, Q$ と置いています。 $H, Q$ はハイパーパラメータで、尤度を最大化するように推定を行います。ハイパーパラメータの推定方法については、野村2016などを参照してください。

このようなモデル設定のもと、観測により得られた値 $Y _ {1:T}=\left \lbrace y_i \right \rbrace _ {i=1}^T$ を用いて未来の観測値の長期予測と状態の推定を行います。

着目するのは以下の事後分布です：

長期予測
- $P(y _ {T+i}|Y _ {1:T})=N(y_{T+i}|\hat{\alpha} _ {T+i}, Z\hat{V} _ {T+i}Z^{\rm T}+H)$
状態推定（平滑化）
- $P(\alpha_t|Y _ {1:T})=N(\alpha_t|\hat{\alpha} _ t, \hat{V} _ t)$

この時の平均ベクトル $\hat{\alpha}$ と共分散行列 $\hat{V}$ は、カルマンフィルタにより逐次的に求めることができます。

ローカルレベルモデル

線型ガウス状態空間時系列モデルの最もシンプルな例として、ローカルレベルモデルを紹介します。ローカルレベルモデルでは、以下のようにブラウン運動する状態 $\alpha_t$ にガウスノイズが加わって観測値 $y_t$ が得られるものと考えます。

${ \begin{aligned} y_t &= \alpha_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, \sigma_\varepsilon^2)\\ \alpha_{t+1} &= \alpha_t + \eta_t, \quad \eta_t \sim N(0,\sigma_\eta^2) \end{aligned} }$

この場合、モデル構造を表すパラメータは $Z=1, T=1, R=1$ のようになります。一方ノイズのパラメータ $H, Q$ は、それぞれノイズの分散 $\sigma _ \varepsilon ^ 2, \sigma _ \eta ^ 2$ に対応します。

上述のローカルレベルモデルで状態の平滑化と長期予測を行うと、以下のような結果が得られます。

f:id:tatamiyatamiyatatatamiya:20210212134144p:plain:w500

上記は馬場2019掲載の架空の売り上げデータ⁴について、観測値の事後分布 $P(y_t|Y _ {1:T})=N(y_t|\hat{\alpha} _ t, \hat{V} _ t+\sigma _ \varepsilon ^ 2)$ を求めたものです。太黒線が観測値、灰線が推定値・予測値 $\hat{\alpha}_t$ です。推定値・予測値周りの影は、推定・予測区間（25~75%, 10~90%）を表しています。

ローカルレベル+季節成分モデル

状態空間モデルでは、1週間・4半期といった周期性のある時系列データを扱うこともできます。

周期性のある時系列は、以下のようにレベル成分( $\mu_t$ )、周期成分( $\gamma_t$ )、観測ノイズ( $\varepsilon_t$ )に分解して扱います。レベル成分というのは、時系列の全体的な変動のことです。

${ \begin{aligned} y_t = \mu_t + \gamma_t + \varepsilon_t \end{aligned} }$

以下では、季節変動を加味した長期予測と、観測値の成分分解を行いたいと思います。そのためには、レベル成分と季節成分を含む状態空間モデルを構築し、次の事後分布を求めます：

季節変動を加味した長期予測
- 観測値: $P(y _ {T+i}|Y _ {1:T})$
観測値の成分分解
- レベル成分（平滑化）: $P(\mu_t|Y _ {1:T})$
- 季節成分（平滑化）: $P(\gamma_t|Y _ {1:T})$

季節成分を含む状態空間モデルを考えます。そのためには、レベル成分 $\mu_t$ ・季節成分 $\gamma_t$ 別々に状態方程式を作ります。今回は、レベル成分 $\mu_t$ は先ほどのローカルレベルモデルで表します。季節成分 $\gamma_t$ については、変動周期を $s$ として以下のように表現します：

${ \begin{aligned} \sum_{i=0}^{s-1} \gamma_{t-i} = \eta^{(\gamma)}_t,\quad \eta^{(\gamma)}_t \sim N(0, \sigma_\gamma^2) \end{aligned} }$

直近 $s$ 個の和が平均的に0であるように置きました。

ここでは周期 $s=4$ とすると、モデルは以下のような形で書けます：

$\alpha_t=(\mu_t, \gamma_t, \gamma _ {t-1}, \gamma _ {t-2}) ^ {\rm T}$ と置くと、上記のモデルを式(\ref{eq:dlm})の形で書き表すことができます。行列 $Z, H, T, R, Q$ はそれぞれ以下のように表されます：

${ \begin{aligned} Z = (1, 1, 0, 0)&, \quad H=\sigma_\varepsilon^2,\\ T= \begin{pmatrix} 1 & -1 & -1 & -1\\ 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \end{pmatrix} , \quad &R= \begin{pmatrix} 1 & 0\\ 0 & 1\\ 0 & 0\\ 0 & 0 \end{pmatrix} , \quad Q = \begin{pmatrix} \sigma_\mu^2 & 0\\ 0 & \sigma_\gamma^2 \end{pmatrix}. \end{aligned} }$

上記のモデルが線型ガウス状態空間モデルに含まれることがわかったので、カルマンフィルタを用いて長期予測・平滑化を行うことができます。

最終的に得られる状態の長期予測・推定分布はガウス分布になります。それぞれ以下のように表すことにします：

${ \begin{aligned} P(\alpha_{T+i}|Y_{1:T})&=N(\alpha_{T+i}|\hat{\alpha}_{T+i}, \hat{V}_{T+i})\\ P(\alpha_t|Y_{1:T})&=N(\alpha_t|\hat{\alpha}_t, \hat{V}_t)\\ \end{aligned}\tag{3}\label{eq:state_smoothed_forecast} }$

これをもとに計算すると、観測値の長期予測値、レベル成分の推定値、季節成分の推定値をそれぞれ以下のように表すことができます：

${ \begin{aligned} P(y_{T+i}|Y_{1:T})&=N(y_{T+i}|Z\hat{\alpha}_{T+i}, Z\hat{V}_{T+i}Z^{\rm T}+\sigma_\varepsilon^2)\\ P(\mu_t|Y_{1:T})&=N(\mu_t|\hat{\alpha}_t^{(1)}, \hat{V}_t^{(1, 1)})\\ P(\gamma_t|Y_{1:T})&=N(\gamma_t|\hat{\alpha}_t^{(2)}, \hat{V}_t^{(2, 2)})\\ \end{aligned}\tag{4}\label{eq:results_forcused_on} }$

$\hat{\alpha}^{(i)}$ は推定状態ベクトル $\hat{\alpha}$ の $i$ 番目の成分、 $\hat{V}^{(i, j)}$ は共分散行列 $\hat{V}$ の $(i, j)$ 番目の成分を表します。

$\hat{\alpha} _ t ^ {(1)}=\hat{\mu}_t, \, \hat{\alpha} _ t ^ {(2)}=\hat{\gamma} _ t$ と置くと、1番目の観測値の状態予測値の平均値は、 $\hat{\mu} _ {T+i} + \hat{\gamma} _ {T+i}$ と表せます。なお、共分散行列の第1項目については、 $Z\hat{V} _ {T+i}Z ^ {\rm T}=\hat{V} _ {T+i} ^ {(1, 1)}+\hat{V} _ {T+i}^{(1, 2)}+\hat{V} _ {T+i} ^ {(2, 1)}+\hat{V} _ {T+i} ^ {(2, 2)}$ と表すことができます。

KFASによるモデル構築

本章ではRを用いて前章で扱ったローカルレベル+季節成分モデルを構築します。

ライブラリはKFASを使います。 KFASは線型状態空間時系列モデル構築に特化したライブラリです。線型ガウスの状態空間モデルを構築し、カルマンフィルタを用いて予測・状態推定を行うことができます。線型ガウス以外にも一部の非ガウスモデルも扱えますが、今回の記事では取り上げません。

今回は、岩波データサイエンスVol.6(伊東2017)の内容をベースに、以下の架空の売り上げデータ⁵を使ってモデリングを行います。

f:id:tatamiyatamiyatatatamiya:20210212134219p:plain:w500

以下では、KFASを用いた基本的なモデルの構築方法と、予測・推定結果の取得方法について解説します。内容の都合上、KFASの利用方法は基本的な部分に留めます。詳細な使い方については参考文献に挙げた書籍や公式ドキュメントを参照してください。

基本的なモデル構築方法

KFASではローカルレベル＋季節成分モデルは、以下のように構築します：

library(KFAS)

# 1. データの読み込み
## 岩波データサイエンスVol.1 松浦(2015)より
## https://github.com/iwanami-datascience/vol1/matsuura/example2/input/data-season.txt
datafile <- "../data/data-season.txt"
data <- read.csv(datafile)
Y <- data$Y

# 2. モデルの定義
## ローカルレベル+季節成分モデル
## H, Qはそれぞれ観測ノイズ、状態ノイズの共分散行列
## H, QにNAを入れることで、要推定パラメータとして扱われる
model <- SSModel(Y ~ SSMtrend(degree = 1, Q=NA) 
                 + SSMseasonal(4, sea.type = "dummy", Q=NA),
                 H=NA)

# 3. ハイパーパラメータ推定
## 2番目の引数は、パラメータ初期値(numeric(3) = (0, 0, 0))
## methodは最適化手法。今回はBFGS(準ニュートン法)。
fit <- fitSSM(model, numeric(3), method="BFGS")

このうち、モデルの定義とハイパーパラメータ推定について簡単に見ていきます。

モデルの定義

# 2. モデルの定義
model <- SSModel(Y ~ SSMtrend(degree = 1, Q=NA) 
                 + SSMseasonal(4, sea.type = "dummy", Q=NA),
                 H=NA)

ここではローカルレベル＋季節成分モデル(\ref{eq:seasonalmodel})を定義しています。ここの部分は観測方程式を以下のように表したものと対応しています：

${ \begin{aligned} y_t \sim N(y_t|\mu_t + \gamma_t, \sigma_\varepsilon^2)\\ \end{aligned} }$

このうちローカルレベル成分 $\mu_t$ の状態方程式はSSMtrend(degree=1, Q=NA)、季節成分 $\gamma_t$ の状態方程式はSSMseasonal(4, sea.type = "dummy", Q=NA)に対応します。観測ノイズの分散 $\sigma _ \varepsilon ^ 2$ に対応する部分はH=NAです。NAが入っている理由は、あとでfitSSMによりハイパーパラメータ推定を行って決めるためです。同様に各成分の状態ノイズについてもQ=NAと置いています。

季節成分を定義したSSMseasonal(4, sea.type = "dummy", Q=NA)について補足します。この第一引数 $4$ は、周期数を指定しています。第二引数にはsea.type = "dummy"とありますが、これによりダミー変数 $\gamma_t$ で季節変動成分を表現するように指定しています。"dummy"のほかに"trigonometric"を指定でき、この場合は三角関数を使って周期性を表現します（詳細は野村2016など）。

ハイパーパラメータ推定

# 3. ハイパーパラメータ推定
fit <- fitSSM(model, numeric(3), method="BFGS")

fitSSMでは、モデルのハイパーパラメータを推定します。今回推定するハイパーパラメータは3つあります: システムノイズの分散 $\sigma _ \varepsilon ^ 2$ 、レベル成分の状態ノイズ分散 $\sigma _ \mu ^ 2$ 、季節成分の状態ノイズ分散 $\sigma _ \eta ^ 2$ です。これらのパラメータについて、尤度が最大になるように最適化を行います。

1番目の引数には、先ほど定義したモデルを指定します。 2番目の引数には、ハイパーパラメータの初期値を指定します。今回はそれぞれ0を入れています。

> numeric(3)
[1] 0 0 0

3番目の引数methodには、最適化手法を指定します。今回はBFGSという、準ニュートン法の一種を用います。 BFGSの意外に使える手法については、最適化ライブラリoptimのドキュメントを参照してください。 KFASはfitSSMの中でoptimを呼んで最適化処理を行っています。

予測・推定結果

推定して得られたハイパーパラメータを用いて、長期予測と状態推定を行います。

長期予測・状態推定ともに、predictを使って行うことができます。 predictは、予測・推定値と予測・推定区間を出力します。予測・推定区間を求めるには、分布がガウス分布となることを思い出しましょう。予測・推定分布が $N(\alpha|\hat{\alpha}, \hat{V})$ のように表せるとします。この時、予測・推定値 $\hat{\alpha}$ の90%区間は分散 $\hat{V}$ を用いて $\hat{\alpha} \pm z(0.05) \sqrt{\hat{V}}$ となります。 $z(a)$ は、標準正規累積分布関数を用いて以下のように定義しました。

${ \begin{aligned} S(x) &= 1-\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x} e^{-t^2/2}dt\\ z(a) &= S^{-1}(a) \end{aligned} }$

区間の範囲はpredictの引数levelで指定することができます。

長期予測

長期予測を行い、

${ \begin{aligned} P(y_{T+i}|Y_{1:T})=N(y_{T+i}|\hat{\mu}_{T+i} + \hat{\gamma}_{T+i}, Z\hat{V}_{T+i}Z^{\rm T}+\sigma_\varepsilon^2) \end{aligned} }$

を求めます。

## 長期予測
predict(fit$model, states='all', n.ahead=8, 
        interval='prediction', level=0.80)

引数について簡単に見ていきます。第一引数にはハイパーパラメータ推定済みのモデルを入れます。第二引数のstatesでは、予測する状態成分を指定しています。今回の場合、レベル成分と季節成分の両方を見るので、allとしています。三番目の引数n.aheadでは、何ステップ先まで予測するかを指定しています。今回は、8期先まで予測します。

四番目の引数のintervalでは、推定区間の計算方法を指定します。 predictionとconfidenceの2種類があり、predictionなら観測ノイズを含む分散 $Z\hat{V} _ {T+i}Z ^ {\rm T}+\sigma _ \varepsilon ^ 2$ をもとに計算を行います。confidenceと指定した場合は観測ノイズを含めず分散 $Z\hat{V} _ {T+i}Z ^ {\rm T}$ をもとに計算します。

今回の対象データに対して8期先までの予測値と50%と80%の予測区間を表したものが、下記のグラフです。

f:id:tatamiyatamiyatatatamiya:20210212134328p:plain:w500

成分ごとの推定

レベル成分

${ \begin{aligned} P(\mu_t|Y_{1:T})&=N(\mu_t|\hat{\mu}_t, \hat{V}_t^{(1, 1)})\\ \end{aligned} }$

を見るには以下のようにします：

# レベル成分
predict(fit$model, interval='confidence', states='level', level=0.80)

推定する状態としてstates='level'を指定
interval='confidence'と指定して観測ノイズを除外する

上記では予測ステップ数n.aheadを指定していません。この場合、入力データに含まれていた全期間を対象に推定が行われます。 n.aheadを指定することで、観測値の長期予測と同様に各レベル成分の長期予測を出力することもできます。

結果は以下のようになります。

f:id:tatamiyatamiyatatatamiya:20210212134303p:plain:w500

同様にして、季節成分

${ \begin{aligned} P(\gamma_t|Y_{1:T})&=N(\gamma_t|\hat{\gamma}_t, \hat{V}_t^{(2, 2)})\\ \end{aligned} }$

を見るためにはstates='seasonal'と指定します。

# 季節成分
predict(fit$model, interval='confidence', states='seasonal', level=0.80)

f:id:tatamiyatamiyatatatamiya:20210212134354p:plain:w500

なお、states='all'もしくはstates=c('level', 'seasonal')と指定すると、観測ノイズを抜いたレベル成分と季節成分の和

${ \begin{aligned} P(\mu_t+\gamma_t|Y_{1:T})&=N(\mu_t+\gamma_t|\hat{\mu}_t + \hat{\gamma}_t, Z\hat{V}_tZ^{\rm T})\\ \end{aligned}\tag{5}\label{eq:components_sum} }$

を取得することもできます。

# レベル成分+季節成分
predict(fit$model, interval='confidence', states='all', level=0.80)

f:id:tatamiyatamiyatatatamiya:20210213124724p:plain:w500

ただし、今回の場合観測ノイズが極めて小さい（ $\sigma ^ 2 _ \varepsilon=0.000183$ ）ため、長期予測のところで得た出力とほとんど変わらない結果が得られます。

# 観測ノイズの分散
fit$model$H
#          [,1]
# [1,] 0.000183

statsmodelsによるモデル構築

以上までの結果を、Pythonを使って再現します。ライブラリはstatsmodelsを用います。

statsmodelsで状態空間時系列モデリングを行うには、UnobservedComponentsというクラスを使います。 UnobservedComponentsクラスを使えば、ローカルレベルモデルや季節成分モデルをはじめとした典型的な状態空間を構築することができます。

以下、KFASの解説と同様に基本的なモデル構築方法、予測・推定結果の取得方法の順で解説します。ただし、状態成分ごとの長期予測と状態成分の和の取得方法については節を分けて解説します。というのも、statsmodelsでこれらの量を取得する方法が、KFASと比較して複雑なためです。

基本的なモデル構築方法

モデルの構築は以下のようにして行います。

import pandas as pd
import statsmodels.api as sm

# 1. データの読み込み
## https://github.com/iwanami-datascience/vol1/matsuura/example2/input/data-season.txt
df_seasonal_sales = pd.read_csv('../data/data-season.txt')
df_seasonal_sales = df_seasonal_sales.rename(columns={'Y': 'sales'})

# 2. モデルの定義
seasonal_model = sm.tsa.UnobservedComponents(df_seasonal_sales['sales'], level='local level', seasonal=4, use_exact_diffuse=True)

# 3. ハイパーパラメータ推定 + 平滑化
seasonal_results = seasonal_model.fit(method='bfgs')

モデルの定義

statsmodelsでは、レベル成分の指定は文字列で行うことができます。ローカルレベルモデルであれば、level='local level'もしくはlevel='llevel'のように指定します。

さらに引数seasonalに周期数を指定することで、ダミー変数による季節成分をモデルに加えることができます。

# 2. モデルの定義
seasonal_model = sm.tsa.UnobservedComponents(df_seasonal_sales['sales'], level='local level', seasonal=4, use_exact_diffuse=True)

引数にuse_exact_diffuse=Trueを入れたは、カルマンフィルタを適用する際の初期値設定方法として散漫初期化(diffuse initialization)を用いるためです。KFASでも散漫初期化を用いているので合わせました。散漫初期化については、野村2016などを参照してください。

なお、ローカルレベルモデル以外のレベル成分の指定方法についてはドキュメントを参照してください。

ハイパーパラメータ推定 + 平滑化

KFASとは異なり、statsmodelsではfitによってハイパーパラメータ推定だけでなく平滑化も行います。

# 3. ハイパーパラメータ推定 + 平滑化
seasonal_results = seasonal_model.fit(method='bfgs')

methodではハイパーパラメータ推定時の最適化手法を指定できます。ここではKFASと同様にBFGSを指定しました。BFGS以外の手法についてはドキュメントのこちらのページを参照してください。

fitを行うと、結果が格納されたオブジェクトが返ってきます。

この結果オブジェクトは、plot_components()というメソッドを持っています。これを使うと、簡単に結果を表示することができます。以下では、

観測データ $Y _ {1:T}$ と一期先予測値 $P(y_t|Y _ {1:t-1})$
レベル成分推定値（平滑化）
季節成分推定値（平滑化）

が表示されます。

fig = seasonal_results.plot_components(which='smoothed')

f:id:tatamiyatamiyatatatamiya:20210212134428p:plain:w500

予測・推定結果

KFASではpredict関数ひとつで長期予測結果も状態推定結果も取得できました。それに対し、statsmodelsでは長期予測と状態推定とで、インターフェースが分かれています。

また、KFASではpredict関数が予測・推定区間まで出力してくれましたが、statsmodelsでは標準偏差を元に自分で計算する必要があります。

長期予測

観測値の長期予測値と分散を得るためには、以下のようにget_forecastメソッドを使います。予測値と分散ともに、返ってきたオブジェクトのプロパティから取得できます。

seasonal_forecast = seasonal_results.get_forecast(steps=8)

seasonal_forecast.predicted_mean
seasonal_forecast.var_pred_mean

KFASと異なり、予測区間を出すには自分で計算してやる必要があります。

from scipy.stats import norm
import numpy as np

critical_value80 = norm.ppf(1. - (1. - 0.80) / 2.)

# upper
seasonal_forecast.var_pred_mean + critical_value80 * np.sqrt(seasonal_forecast.var_pred_mean)

# lower
seasonal_forecast.var_pred_mean - critical_value80 * np.sqrt(seasonal_forecast.var_pred_mean)

こうして得られた結果をプロットしたものが以下になります：

f:id:tatamiyatamiyatatatamiya:20210212134453p:plain:w500

成分ごとの推定

平滑化により得られたレベル成分と季節成分の推定値( $\hat{\mu} _ t, \hat{\gamma} _ t$ )・分散( $\hat{V} _ t ^ {(1,1)}, \hat{V} _ t ^ {(2,2)}$ )は、fitで返ってきたオブジェクトのプロパティから取得できます。

# レベル成分
seasonal_results.level['smoothed']
seasonal_results.level['smoothed_cov']

# 季節成分
seasonal_results.seasonal['smoothed']
seasonal_results.seasonal['smoothed_cov']

こちらも推定区間の算出は自分で行う必要があります。

KFASと比べて取得しにくい量

以上のように、statsmodelsでも簡単に状態空間モデリングを行うことができました。結果の取得もKFASにやや劣るものの比較的容易でした。

しかし、KFASではpredict関数ひとつで簡単に取得できた量を、statsmodelsでは容易に見つけられないケースもあります。

例えば、下記の二点があります:

成分ごとの長期予測値と分散
レベル成分と季節成分の和の予測・推定値と分散

これらをstatsmodelsで取得するには、初見ではわからない場所に格納されているオブジェクトを探しに行かなくてはいけません。さらにその際、予測・推定結果がベクトル・行列でどのように表現されるかを知っておかなくてはいけません。

以下では、上記2種類の量の取得方法を解説したうえで、KFASで出力したものと同じ図を作成します。

成分ごとの長期予測

KFASでは各成分ごとの長期予測値を、観測値の長期予測と同様にpredictを使って取得することができました。しかし、statsmodelsではそもそも成分ごとの長期予測を得るための関数・メソッド、プロパティが存在しません。そのため、式(\ref{eq:state_smoothed_forecast})で解説した状態ベクトル $\hat{\alpha} _ {T+i}$ とその共分散行列 $\hat{V} _ {T+i}$ を直接探しにいきます。

予測時の生の計算結果を取得するには、get_forecast()で返ってきた予測結果オブジェクトforecastの持つforecast.prediction_results.resultsという内部オブジェクトを見ます。この内部オブジェクトの持つプロパティpredicted_stateとpredicted_state_covが、それぞれ状態ベクトル $\hat{\alpha} _ {T+i}$ と共分散行列 $\hat{V} _ {T+i}$ に相当します。

predicted_stateとpredicted_state_covからレベル成分の予測値と分散を取得します。レベル成分の予測値 $\hat{\mu} _ {T+i}$ は、状態ベクトル $\hat{\alpha} _ {T+i}$ の第一成分でした（式(\ref{eq:results_forcused_on})参照）。したがって、predicted_stateの0番目の成分がレベル成分の予測値 $\hat{\mu} _ {T+i}$ です。同様にしてレベル成分の分散 $\hat{V} _ {T+i} ^ {(1,1)}$ は、predicted_state_covの $(0, 0)$ 成分です。

# レベル成分長期予測
## 平均値
seasonal_forecast.prediction_results.results.predicted_state[0]

## 分散
seasonal_forecast.prediction_results.results.predicted_state_cov[0][0]

先ほどの平滑化推定量と合わせてプロットすると、以下のようになります：

f:id:tatamiyatamiyatatatamiya:20210212134522p:plain:w500

同様にして季節成分の予測値 $\hat{\gamma} _ {T+i}=\hat{\alpha} _ {T+i} ^ {(2)}$ と分散 $\hat{V} _ {T+i} ^ {(2,2)}$ も以下のように取得できます。

# 季節成分長期予測
## 平均値
seasonal_forecast.prediction_results.results.predicted_state[1]

## 分散
seasonal_forecast.prediction_results.results.predicted_state_cov[1][1]

f:id:tatamiyatamiyatatatamiya:20210212134543p:plain:w500

成分和の予測・推定

KFASではpredict関数の引数でstates=c('level', 'seasonal')とすることで、レベル成分と季節成分両方の和について予測・推定値と分散を出力できました。このような任意の成分を組み合わせた予測・推定をstatsmodelsで行うにも、予測値・推定値の共分散行列を探しにいく必要があります。以下では、レベル成分と季節成分の和について平滑化推定と長期予測の結果を取得します。

まずはレベル成分+季節成分の平滑化推定値と分散を求めます。式(\ref{eq:components_sum})から、推定値は $\hat{\mu} _ t + \hat{\gamma} _ t$ 、分散は $Z\hat{V} _ tZ ^ {\rm T}=\hat{V} _ t ^ {1, 1}+\hat{V} _ t ^ {1, 2}+\hat{V} _ t ^ {2, 1}+\hat{V} _ t ^ {2, 2}$ と表せます。

推定値については個々の成分推定値の和をとるだけです。 $\hat{\mu} _ t, \hat{\gamma} _ t$ はレベル成分と季節成分の平滑化推定値ですから、前の節のようにfit結果オブジェクトから直接取得できます。

# レベル成分+季節成分 平滑化推定値
seasonal_results.level['smoothed'] + seasonal_results.seasonal['smoothed']

一方で分散については、推定値の共分散行列 $\hat{V} _ t$ を探さなくてはいけません。そのためには、fit結果オブジェクトもつfilter_resultsを見ます。このオブジェクトの持つsmoothed_state_covというプロパティが共分散行列 $\hat{V}_t$ です。

# レベル成分+季節成分 平滑化分散

## NG:
## seasonal_results.level['smoothed_cov'] + seasonal_results.seasonal['smoothed_cov']

seasonal_results.filter_results.smoothed_state_cov[:2, :2].sum(0).sum(0)

同様に、レベル成分+季節成分の長期予測結果を取得します。こちらについては、成分ごとの長期予測と同様にget_forecastで返ってきたオブジェクトから予測状態ベクトル $\hat{\alpha} _ {T+i}$ と共分散行列を $\hat{V} _ {T+i}$ を取得して計算を行います。

# レベル成分+季節成分 長期予測

## 予測値
seasonal_forecast.prediction_results.results.predicted_state[0] + seasonal_forecast.prediction_results.results.predicted_state[1]

## 分散
seasonal_forecast.prediction_results.results.predicted_state_cov[:2, :2].sum(0).sum(0)

以上を合わせてプロットすると、以下のようになります：

f:id:tatamiyatamiyatatatamiya:20210212134608p:plain:w500

まとめ

本記事ではPythonで状態空間時系列モデリングを行う方法を理解するために、ローカルレベル+季節成分モデルの構築をR→Pythonの順で行いました。

PythonでもstatsmodelsのUnobservedComponentsクラスを使うことでRのライブラリKFASと同様のモデリング結果を出力することができました。ただし、statsmodelsでは予測・推定量の取得方法が分かりにくいという難点がありました。具体的には、

予測・推定結果を取得するためのインターフェイスが統一されていない
求めたい量によっては、状態ベクトルや共分散行列から直接計算を行わないといけない

ということが挙げられます。予測・推定量の取得方法をまとめると、以下の表のようになります：

	statsmodels	KFAS（予測・推定値±予測・推定区間）
状態平滑化 $P(\alpha_t \mid Y _ {1:T})=N(\alpha_t \mid \hat{\alpha} _ t, \hat{V} _ t)$	$\hat{\alpha} _ t$ : `UnobservedComponentsResults.filter_results.smoothed_state` / `UnobservedComponentsResults.{level/seasonal}['smoothed']` $\hat{V} _ t$ : `UnobservedComponentsResults.filter_results.smoothed_state_cov` / `UnobservedComponentsResults.{level/seasonal}['smoothed_cov']`	`predict(fit$model, states='{all/level/seasonal}', interval='confidence')`
状態長期予測 $P(\alpha _ {T+i} \mid Y _ {1:T})=N(\alpha _ {T+i} \mid \hat{\alpha} _ {T+i}, \hat{V} _ {T+i})$	$\hat{\alpha} _ {T+i}$ : `PredictionResults.prediction_results.results.predicted_state` $\hat{V} _ {T+i}$ : `PredictionResults.prediction_results.results.predicted_state_cov`	`predict(fit$model, states='{all/level/seasonal}', interval='confidence', n.ahead=10)`
観測値平滑化 $P(y_t\mid Y _ {1:T})=N(y_t\mid Z\hat{\alpha} _ t, Z\hat{V} _ t Z ^ {\rm T}+\sigma _ \varepsilon ^ 2)$	推定値: `UnobservedComponentsResults.filter_results.smoothed_forecasts` 分散: `UnobservedComponentsResults.filter_results.smoothed_forecasts_error_cov`	`predict(fit$model, states='{all/level/seasonal}', interval='prediction')`
観測値長期予測 $P(y_{T+i}\mid Y _ {1:T})=N(y _ {T+i}\mid Z\hat{\alpha} _ {T+i}, Z\hat{V} _ {T+i} Z ^ {\rm T}+\sigma _ \varepsilon ^ 2)$	予測値: `PredictionResults.predicted_mean` 分散: `PredictionResults.var_pred_mean`	`predict(fit$model, states='{all/level/seasonal}', interval='prediction', n.ahead=10)`

KFASではpredict関数一つで予測・推定区間の取得まで行うことができます。現状のstatsmodelsは使いづらい部分が目立ちますが、インターフェイスを改良してドキュメントを整えればそれなりに使いやすくなるのではないかな、と個人的には考えています。

参考資料

本記事の図を作成する際に用いたR, Pythonのコードは、こちらにアップしています。

畳庵〜tatamiya practice〜

機械学習・統計モデリング練習帳。読書・source code読みの記録や，実装など。

statsmodelsによる状態空間モデリングを、KFASから「翻訳」して学ぶ（R経由で学ぶPythonの状態空間時系列モデリング#2）

はじめに

忙しい人のために

状態空間モデルによる季節変動のある時系列データの扱い方

線型ガウス状態空間モデル（前回の復習）

ローカルレベルモデル

ローカルレベル+季節成分モデル

KFASによるモデル構築

基本的なモデル構築方法

モデルの定義

ハイパーパラメータ推定

予測・推定結果

長期予測

成分ごとの推定

statsmodelsによるモデル構築

基本的なモデル構築方法

モデルの定義

ハイパーパラメータ推定 + 平滑化

予測・推定結果

長期予測

成分ごとの推定

KFASと比べて取得しにくい量

成分ごとの長期予測

成分和の予測・推定

まとめ

参考資料