渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

マルコフ状態転換モデルのRパッケージ{MSwM}の使い方（異常値検出・ステータス変化検出などに有用）

統計学 R 時系列分析

CRANパッケージ{MSwM}の大体の使い方が分かったので簡単に共有します。


なお、しつこいようですがマルコフ状態転換モデルについてはこのブログではすっかりお馴染みの以下のテキストをご参照のこと*1。僕もまだ勉強中です。


経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

Time Series Analysis

Time Series Analysis


なお、沖本テキストでは簡単な説明のみに留められているので、Hamiltonテキストの方がお薦めです。ただし結構難解なので要注意。


マルコフ状態転換モデル(Markov switching model / Regime-switching model)について


マルコフ状態転換モデルの原理的説明については、ここでは大幅に割愛します*2。基本的なアイデアを簡潔に書いておくと、

  1. 観測時系列がn個の状態に分かれ得ると仮定して、パラメータフリーなAR過程をn個置く
  2. 観測時系列はn \times nの遷移確率行列で表されるマルコフ過程に従って状態変化すると仮定する
  3. EMアルゴリズムを用いて、観測時系列がn通りある状態のいずれかであったかの事後確率列を求める


というものです。要は、AR過程に対する隠れマルコフモデルと思ってください。EMアルゴリズム周りについては、最新の理論を当てはめれば今後はもっと良くなるのかもしれません。

具体的な{MSwM}パッケージの使い方


とりあえずVignetteを読んでみて下さい。{MSwM}はマルコフ状態転換モデルを推定するに当たって、何かしらの説明変数を仮定してそれとの回帰モデルに基づいて、EMアルゴリズムで事後確率を求めるというやり方をしています。なので、いかなる時系列であれ、（どんなに恣意的でも良いので）説明変数を持ってくる必要があります。ここは注意が必要。手順としては、

  1. 事後確率列を求めたい時系列yを持ってくる
  2. 何でも良いので説明変数になり得る時系列x_1, x_2, \cdots, x_nを持ってくる
  3. yとx_nとでデータフレームを作り、lm(){stats}関数で回帰モデルmodelを推定する
  4. msmFit(){MSwM}関数にmodel、状態数k、推定パラメータリストswを与えて、マルコフ状態転換モデルを推定する


という感じの流れになります。クセモノなのがx_nで、特に何も思い付かない時は単なる単調増加関数（ってかただの直線）を与えても良いのですが、y次第ではある程度何がしかの「参考」になりそうなデータを与えてやった方がベターです。


もっと平たく言えば、事後確率列を求めたい時系列（スペンドUUなど）に対して、明らかに連動して変化していると思しき時系列（DAUや新規UUなど）が別にある場合は、それをx_nとして与えた方が確実ということです。

異常値検出のケース


まずサンプルデータを作ります。

> x1<-matrix(rnorm(300)) # ホワイトノイズ
> x2<-matrix(rep(c(15,rep(0,49)),6)) # スパイク異常値が6回生じる時系列
> x<-x1+x2 # これらを足し合わせる
> xdum<-matrix(1:300) # 説明変数となる、直線状のダミー変数を作る
> xmsm<-data.frame(cbind(x,xdum)) # データフレームにまとめる
> names(xmsm)<-c("datx","dumx") # 名前をつける


f:id:TJO:20130620154222p:plain


こんな感じのデータが出来上がります。こいつに対して、以下のような手順でマルコフ状態転換モデルを推定し、異常値を検出してみます。

> modelx<-lm(datx~.-1,xmsm) # 回帰モデルを求める
> msmModelx<-msmFit(modelx,k=2,sw=rep(T,2))
# マルコフ状態転換モデルを推定する。kは状態数、
# swは推定したいパラメータ数をlogicalのT / Fで指定する
> plotProb(msmModelx,which=1) # 全体の事後確率列のリストをまとめてプロット
> plotProb(msmModelx,which=2) # 1番目の事後確率列をプロット
> plotProb(msmModelx,which=3) # 2番目の事後確率列をプロット


まず最初に、全体の事後確率列のリストをプロットしてみます。異常値とそうでないところが分離されているらしい、というのが分かります。

f:id:TJO:20130620154315p:plain

次に、事後確率列のプロットを原時系列プロットの上に陰影をつける形で表示してみます。

f:id:TJO:20130620154326p:plain

f:id:TJO:20130620154334p:plain


きちんとスパイク異常値がマルコフ状態転換モデルによって、通常値から分離してモデリングされていることが分かりますね。


ステータス変化検出のケース


これもサンプルデータ作りから始めます。

> y1<-matrix(rnorm(150)-5) # ベースラインが低いホワイトノイズ
> y2<-matrix(rnorm(150)+5) # ベースラインが高いホワイトノイズ
> y<-rbind(y1,y2) # 2つを合わせる
> dumy<-rbind(matrix(1:150),matrix(150:1)) # ダミー変数として「山」の字を持ってくる
> ymsm<-data.frame(cbind(y,dumy)) # データフレームにする
> names(ymsm)<-c("daty","dumy") # 名前をつける


f:id:TJO:20130620154920p:plain


これまたこんな感じのデータが出来上がります。「どう見ても説明変数が恣意的じゃねーかゴルァ」とか言わないように（笑）。ともあれ、このデータに対してマルコフ状態転換モデルを推定してみます。

> modely<-lm(daty~.-1,ymsm)
> msmModely<-msmFit(modely,k=2,sw=rep(T,2))
> plotProb(msmModely,which=1) 
> plotProb(msmModely,which=2)
> plotProb(msmModely,which=3)


全体の事後確率列のリストを見ると、綺麗に「カックン」とベースラインが変わったところが検出されているのが分かりますね。

f:id:TJO:20130620155114p:plain

陰影プロットで見れば、もちろん一目瞭然です。

f:id:TJO:20130620155213p:plain

f:id:TJO:20130620155220p:plain


同様に、途中でベースライン変動というステータス変化があった場合でも、マルコフ状態転換モデルでバッチリ検出できるということが分かりました。


ちなみに、この例で異常値検出の時と同じように直線状の説明変数xを与えると、ものの見事に失敗します*3。思った以上にx_nの設定が大事っぽいです。


補足


もちろん観測時系列yが、「異常値」と「ステータス変化」の両方を同時に含むケースもあります。この場合は、連動していると思しき別の時系列を見繕ってx_nとして置いてやって、k \geq 3とすればきちんと異常値とステータス変化を同時にモデリングすることが可能です。


ただ、それができるかどうか？というのは今のところ完全に経験（暗黙知）に依存している部分が大きいです。この辺はまた改めてマルコフ状態転換モデルについて書くときに、きちんと論じましょうということで。。。

*1:1万回でも10万回でも書きますが、アフィリエイトの類はやってませんよー

*2:また改めて「Rで計量時系列分析」シリーズと称して詳しく書く予定です

*3:何故か「ずーっと平坦」という推定結果になる