この記事はStan Advent Calendar 2018

の23日目の記事です。本記事では状態空間モデルを用いた時系列予測について実データでの分析例を紹介します。

2015年度人工知能学会全国大会（第29回）での論文、

状態空間モデルを用いた検索トレンドとページビューからの自動車販売台数の予測, 角田孝昭, 吉田光男, 津川翔, 山本幹雄

www.jstage.jst.go.jp

と同内容の状態空間モデルによる予測をStanを用いて行ってみました。この論文では、自動車の販売台数の月次データの予測を状態空間モデルを用いて行っていて、説明変数として絵googleの検索量を用いることで予測が改善するかどうかを調べています。

使用データについて

今回使用したデータは論文と同様に日本自動車販売協会連合会

自販連のホームページ

からとってきた車種ごとの月次販売台数データと、Googleトレンドからとってきた検索量のデータです。販売台数データについてはpdfファイルから

suryu.me

を参考に取り出しました。データを取り出す過程もいろいろ面倒なところがあったのですが、これについては記事を改めて紹介したいと思います。

予測モデルについて

状態空間モデルについては本Advent Calenderでの

logics-of-blue.com

をご覧ください。簡単に紹介すると、時系列のようなデータに対して、その観測値を直接予測するのではなく状態という観測できない隠れた変数を用意し、状態についての時系列モデルに加えて状態にノイズを乗せて観測値が得られるというモデルです。

もっともシンプルなモデルとしては、

$\displaystyle{ y_t = \mu_t + v_t\\ \mu_t = \mu_{t-1} + w_t }$

というものがあります。

ここで、

tが時刻
y_tが観測値
v_tが観測誤差
μ_tが状態
w_tが状態の変動

を表します。このw_tやv_tが正規分布に従うと仮定すると、μはランダムウォークになります。つまり、一時点前の状態μ_t-1に正規分布から発生した量w_tが加わることで、時刻tでの状態μ_tが得られるというのが状態のモデルです。さらにこの状態から、正規分布から発生した誤差が加わったことで観測値y_tが得られるというのが状態空間モデルです。推定するパラメータは状態μ_t、及び観測誤差w_tや状態変動v_tの分散で、これらを推定したうえでy_tを予測します。

今回のコードについては上記の松浦さんの本に同じような状態空間の例がありましたので、そちらを参考にさせていただきました。本で用いられているコードは

github.com

にありますので、そちらもご参照ください。また結果の図示については

statmodeling.hatenablog.com

を参考にしました。（もし問題がある場合にはお知らせ下さい。）

元論文ではいくつかの車種について実験を行なっていますが、今回は販売台数一位のノートのみを用いています。

まずは論文のベースラインモデルですが、2階差分トレンドで12ヶ月周期の季節項入りのモデルを、松浦さんの本にちょうど同じものがあったので、参考にしました。

モデルの式は以下の通りです。

$\displaystyle{ y_t = \mu_{y,t} + \gamma_{t} + v_t\\ \mu_{y,t}=2\mu_{y,t-1}-\mu_{y,t-2}+w_t }$

で、

y_tが販売台数の観測値
γ_y,tが販売台数の季節成分
μ_y,tが販売台数のトレンド
v_tが販売台数の観測誤差
w_tが販売台数のトレンドの変動

です。

次に提案モデルですが、こちらは検索量も2階差分の季節成分入りのトレンドモデルを用いて状態を予測し、販売量の状態に一次式で説明変数として検索量のトレンド成分を加えるというものです。

モデルは少々複雑ですが以下の通りです。

$\displaystyle{ y_t = \mu_{y,t} + \gamma_{y,t} + b \mu_{x,t}+v_{y,t}\\ \mu_{y,t}=2\mu_{y,t-1}-\mu_{y,t-2}+w_{y,t}\\ x_t = \mu_{x,t} + \gamma_{x,t} + v_{x,t}\\ \mu_{x,t}=2\mu_{x,t-1}-\mu_{x,t-2}+w_{x,t} }$

これは検索量についての状態空間モデルと販売台数についての状態空間モデルの組み合わせです。

x_tが検索量の観測値
γ_x,tが検索量の季節成分
μ_x,tが検索量のトレンド
v_x,tが検索量の観測誤差
w_x,tが検索量のトレンドの変動

というのが下二つの式に現れる変数で、こちらが検索量についての季節成分入り2階差分トレンドモデルになっています。

y_tが販売台数の観測値
γ_y,tが販売台数の季節成分
μ_y,tが販売台数のトレンド
bが販売台数への検索量の影響の重み
v_y,tが販売台数の観測誤差
w_y,tが販売台数のトレンドの変動

というのが上二つの式に現れる変数です。こちらも同様に季節成分入り2階差分トレンドモデルですが、販売台数の観測値の式に検索量の状態が説明変数として加わっていることに注意してください。

モデル部分のStanコードはこんな感じです

model {
  alpha_X[3:N] ~ normal(2*alpha_X[2:(N-1)] - alpha_X[1:(N-2)], s_a_X); //検索量トレンドの状態モデル
  for(t in 12:N){
    season_X[t] ~ normal(-sum(season_X[(t-11):(t-1)]), s_season_X); //検索量季節成分
  }
  X ~ normal(x_mean, s_X); //検索量の観測モデル
  alpha_Y[3:N] ~ normal(2*alpha_Y[2:(N-1)] - alpha_Y[1:(N-2)], s_a_Y); //販売台数トレンドの状態モデル
  for(t in 12:N){
    season_Y[t] ~ normal(-sum(season_Y[(t-11):(t-1)]), s_season_Y); //販売台数季節成分
  }
  Y ~ normal(y_mean + b * alpha_X, s_Y); //販売台数の観測モデル
}

またそのバリエーションとして、検索量を同じ月のものではなく1月前にしたもの、2期前にしたもの、また2期前までの3期分を使ったもの、の3種類のモデルを考えます。

モデルの式はそれぞれ以下の通りです。

$\displaystyle{ y_t = \mu_{y,t} + \gamma_{y,t} + b \mu_{x,t-1}+v_{y,t}\\ \mu_{y,t}=2\mu_{y,t-1}-\mu_{y,t-2}+w_{y,t}\\ x_t = \mu_{x,t} + \gamma_{x,t} + v_{x,t}\\ \mu_{x,t}=2\mu_{x,t-1}-\mu_{x,t-2}+w_{x,t} }$

このモデルでは、一期前の検索量のトレンド成分を説明に用いるので、μ_x,t-1がy_tの式に入っています。

$\displaystyle{ y_t = \mu_{y,t} + \gamma_{y,t} + b \mu_{x,t-2}+v_{y,t}\\ \mu_{y,t}=2\mu_{y,t-1}-\mu_{y,t-2}+w_{y,t}\\ x_t = \mu_{x,t} + \gamma_{x,t} + v_{x,t}\\ \mu_{x,t}=2\mu_{x,t-1}-\mu_{x,t-2}+w_{x,t} }$

このモデルでは、二期前の検索量のトレンド成分を説明に用いるので、μ_x,t-2がy_tの式に入っています。

$\displaystyle{ y_t = \mu_{y,t} + \gamma_{y,t} + b_0 \mu_{x,t} + b_1 \mu_{x,t-1} + b_2\mu_{x,t-2}+v_{y,t}\\ \mu_{y,t}=2\mu_{y,t-1}-\mu_{y,t-2}+w_{y,t}\\ x_t = \mu_{x,t} + \gamma_{x,t} + v_{x,t}\\ \mu_{x,t}=2\mu_{x,t-1}-\mu_{x,t-2}+w_{x,t} }$

このモデルでは、二期前までの検索量のトレンド成分を説明に用いるので、μ_x,t, μ_x,t-1, μ_x,t-2の三つがy_tの式に入っています。

予測結果について

学習データとして2014/1から2017/6までのデータを使い、2017/7から2018/6までをテストデータとしています。これを用いてStanによりMCMCサンプリングした結果のRhatをmcmc_rhatでプロットしたもの以下の通りです。

f:id:unaoya:20181223235814p:plain — Rhat

f:id:unaoya:20181223232429p:plain — Rhat

複数期のトレンド成分を用いたモデルは収束が悪いですが、今回はここを改善する時間がなかったのでこのままとします。

また時系列での予測のプロットは以下の通りです。

f:id:unaoya:20181223235811p:plain

f:id:unaoya:20181223235834p:plain

f:id:unaoya:20181223235840p:plain

f:id:unaoya:20181223235850p:plain

f:id:unaoya:20181223235822p:plain

太い線が実際のデータで細い線がサンプリングした結果の中央値です。また, 色の濃い部分が10%から90%の範囲、薄い部分が2.5%から97.5%です。ベースラインモデル以外では、学習データの部分で実際の値と予測値の中央値が少しずれている部分があるのが気になります。

これらのモデルの予測の評価を行いましょう。元論文ではカルマンフィルタにより予測値を計算しそれと実際の値の誤差二乗和を計算していますが、今回はStanでサンプリングした予測値の平均値とテストデータとの誤差二乗和で評価を行います。また、元論文では相対的な誤差を計算していますが、今回は1車種しか予測しないので、そのままの値で計算します。

それぞれのモデルで誤差を計算した結果は以下の通りです。

ベースラインモデルが349,513,623
提案モデルが306,010,351
一期前モデルが125,867,523
二期前モデルが397,689,428
複数期モデルが1,633,516,793

となっており、一期前の検索量を入れたモデルがもっともよく予測できているという結果でした。

他の車種も含めた分析であったり、他のモデルを試してみたり、またモデル比較についても交差検証やWAIC, WBICその他の指標を用いるなどいくつか試してみたいことはありますので、また実験してみたいと思います。最後までお付き合いありがとうございました。

なお、今回分析で用いたコードは以下にあります。

検索量を用いた状態空間モデルによる売上予測

もし誤りや、改善点などお気付きでしたらお知らせいただけると幸いです。

sigma

統計に関すること

検索量を用いた状態空間モデルによる売上予測

使用データについて

予測モデルについて

予測結果について