Home＞Ｒな 予測＞BUGSで状態空間モデル

　季節とトレンド

★新サイト完成しました!
３秒後に自動的に移動します

変わらない方は こちらからどうぞ
http://logics-of-blue.com/%E5%AD%A3%E7%AF%80%E3%81%A8%E3%83%88%E3%83%AC%E3%83%B3%E3%83%89/

　dlmパッケージを使って、ちょっと複雑な正規線形状態空間モデル（動的線形モデル）を推定します。
　今回は、トレンド成分と季節成分の入った状態空間モデルを推定し、予測を行うところまでを行います。

状態空間モデル関連のページ
　 状態空間モデル　　　　　　　状態空間モデルのことはじめ　
   ｄｌｍの使い方 　　　 　　　　　Ｒで正規線形状態空間モデルを当てはめる
　　ロー カルレベルモデル　　　ｄｌｍパッケージを使ってローカルレベルモデルを当てはめる
　　季 節とトレンド　　　　　　　　dlmパッケージを使って季節成分とトレンドの入ったモデルを作る

　　　　　　　　　　　　　　　　　　　　　　 　　　　　　　　　　　　　　　作成日　　２０１２年 ８月 ２１日　
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　最終更新　２０１２年８月２１日

ローカルレベルモデルを発展させる

　こちらで書いたように、一番単純な動的線形モデルであるロー カルレベルモデルは、データから正しそうな状態を推定するということはある程度可能でしたが、予測には何の役にも立ちませんでした。
　これは、モデルを作るのに使った「データに対する知識」があまりにも少なかったことが原因です。

　たとえば、「このデータは上昇（下降）トレンドがある」とか、「毎年８月には高い値を示す」といった季節変動なんかの知識を取り込むとより良い予 測が可能です。

　というわけで、今回は、その　「トレンド」　と　「季節」　に焦点を当てて、これらの要因を動的線形モデルに付け加える方法について記します。

トレンド

トレンドは大きく分けて２つあります。
　一つは、確定的トレンド
　もう一つが確率的トレンド
です。

確定的トレンドは、回帰分析を思い出してもらえるとよいです。散布図に 直線を引っ張るアレです。
横軸に時間を取って、縦軸にたとえばビールの売り上げを置いたら、それこそ右肩上がりに一定のスピードで増加していく、というのが確定的トレンド。
　図で示すとこんな感じです。

　上のグラフの場合は、時間が１増えるごとに、ビールの売り上げは２ずつ増えていっていることがわかると思います。
　一定して、時間が１増えると、売り上げが２増える、というのが、確定的トレンドです。

一方、確率的トレンドは、その名の通り、確率的に増加(減少)分が変化します。
この変化のパターンもランダムウォークです。
一期前のトレンドと似たようなトレンドになるけれども、ちょっと値が変わるということです。

ローカル線形トレンドモデル

　トレンドとローカルレベルモデルを合体させたもののイメージをつかんでいただくために、例を使って説明します。

２０００年におけるビールの売り上げは１０本。その時のトレンドは、＋４でした。
順当にいけば２００１年は１４本売れることになります。

　最尤法をした結果、観測誤差は±７の範囲を持つことがわかり、状態の変動の大きさは±４であることが推定できたとします。さらに、トレンドの変動 の大きさは±２でした。

　２０００年の状態は１０でした。
　２００１年は、１０から±４の範囲内で動く可能性があります。さらに、トレンドは４±２なので、２〜６の範囲で変動します。
　サイコロを振ったら状態の変化は＋３になったとします。さらに、トレンドの変化は−２になったとします。
　すると、２００１年の状態は　

１０（元の値）　＋　３（状態の変動）　＋　４−２（トレンド）　＝　１５　　　となります。

　ノイズの大きさは常に±７なので、２００１年の観測値は、８〜２２の範囲内に収まるだろうと考えられます。　

　予測を行う際、サイコロは正負どちらに転ぶかわかりません。そういう時は「変化しない」を前提として予測するんでした。ですから、２００１年の予 測値は

１０（２０００年の状態の値）　＋　０　＋　４（２０００年のトレンドの値）　＝　１４

という風に予測されます。で、２００１年のデータが手に入ったら、それをもとにして予測された状態を補正する、という流れは一緒です。

季節成分

　今度は季節成分について説明します。
　季節成分の入れ方には三角関数（サインとかコサインとか）を使ったやり方と、ダミー変数を使ったやり方の2種類がありますが、ここでは後者のダミー変数 を使う方法だけを説明します（dlmパッケージを使って三角関数な季節成分を推定することも可能です）。

　ダミー変数を使ったやり方とは、単純にひと月ごとに売り上げの大きさの補正をすることです。1月にはビールは売れないけど、8月にはよく売れる、 といったような感じ。
　この季節成分にもシステムノイズを入れ込むことは可能です。ノイズがあるのが確率的季節成分、ないのが確定的季節です。
　たとえば、1月には売れ行きが悪くて−５の補正、2月には若干ましになって−３の補正を受けたとします。
　２０００年1月の売り上げは１０でした。季節成分だけを見ると、２月は1月よりプラス２になっているので、２０００年２月の売り上げは１２となるはずで す。これがシステムノイズがない場合。
　システムノイズがあるときは、ローカルレベルモデルなんかと同じく、季節ごとにさいころを投げてやります。季節には±１の変動があったとしたら、 ２０００年２月の値は、１１〜１３ということになります。
　もちろんこの季節のブレだけでなく、ローカルレベルモデルで表された状態の変動やトレンドの変動なんかも加えていくことが可能です。

ｄｌｍで季節とトレンド

今回使ったソースはこちらにおいてあります。　

今回は、季節成分とトレンド両方がありそうなデータとして、AirPassengersというデータを使います。
　もともとのデータを対数変換したものを解析の対象として使います。

> data <- log(AirPassengers)
> plot(data,type="o")

　まずは、状態の変化に一切確率的な変動がない、プロセスノイズが全部０なモデルを推定します。
　ローカルレベルモデルで示したよう、システムノイズが０な状態空間 モデルは、普通の線形回帰をした場合と一緒の結果が出ました。ですので、このモデルは重回帰分析を使うことによっても推定可能です。暇な方は試してみてく ださい。

ｄｌｍパッケージを落としてから

library(dlm)

あとは、以前に示したように４ステップで推定完成です。

#　Step1
#　モデル作成のための関数を作る
build.4<-function(theta){
    dlmModPoly(order=2,dV=exp(theta[1]),dW=c(0,0))+
    dlmModSeas(fr=12,dW=c(0,rep(0,10)),dV=0)
    }

#　Step2
#　MLEでパラメタ推定。
fit.4<-dlmMLE(data,parm=dlmMLE(data,parm=c(1),build.4,method="Nelder-Mead")$par,build.4,method="BFGS")

# 推定されたパラメタを使ってモデルを作り直す
DLM.4<-build.4(fit.4$par)

# Step3
# カルマンフィルター
Filt.4<-dlmFilter(data,DLM.4)

# Step4
# スムージング
Smooth.4<-dlmSmooth(Filt.4)

モデルの型を決める作業が一番大きく変わったと思います。

build.4<-function(theta){
    dlmModPoly(order=2,dV=exp(theta [1]),dW=c(0,0))　+
    dlmModSeas(fr=12,dW=c(0,rep(0,10)),dV=0)
    }

において、緑の部分がローカル線形トレンドモデル（ローカルレベルモデル＋トレンド）で、紫が季節の成分です。両者を　＋　でつないでいることがわ かると思います。

　緑の部分から説明します。黒い太線になっているorder=2ですが、これが１ならば単なるローカルレベルモデ ルで、２にするとトレンド付きになります。ちなみにデフォルトは２なので、いちいち書かなくても実は問題ありません。
　ｄVは観測誤差で、この大きさは推定しています。でも、状態そのものは確率的に変動しないというモデルを今回は組んでいるので、dWは０を突っ込んでい ます。ここで、dWに二つも０を入れていますが、一つ目は（ローカルレベルモデルで説明した）状態の変動の大きさを表し、二つ目はトレンドの変動の大きさ を指定しています。ですので、片方だけ０とか言ったモデルももちろん推定できます。

　紫の部分は季節を表したところですが、これにはdlmModSeas()を使います。fr＝１２とは、一年の間に１２こデータがありますよという 指示。たとえば春夏秋冬４つの区分しかなければｆｒ＝４とします。
　ｄWにやけにたくさん０があります（最初の０に追加してrep()でさらに０を１０個生産している）が、これは最後の１０個の０はほぼ確定です。気にせ ず０を入れればよいです。
　ｄWの一番最初は、季節が確率的変動をする場合は、この大きさも推定することになります。
　観測誤差のほうはdlmModPoly()で推定されているので、ここには０を突っ込んでおきます。

　推定された観測誤差は

> DLM.4$V
          [,1]
[1,] 0.0035172

となりました。ふつうにDLM.4と推定されたモデルの型を表示させるとすご くたくさんの数値が出てきますが、基本的な見方はローカルレベルモデルと一緒なので省略します。季節があるから、その分表示する量が多くなったんですね。

　グラフを書きます

plot(data,col=1,type="o")
lines(dropFirst(Filt.4$m)[, 1] + dropFirst(Filt.4$m)[, 3],col=2)
lines(dropFirst(Smooth.4$s)[, 1] + dropFirst(Smooth.4$s)[, 3], col=4)
legend("bottomright", pch=c(1,NA,NA),
    col=c(1,2,4), lwd=1, legend=c("data","Filter","Smooth"))

　lines(dropFirst(Filt.4$m)[, 1] + dropFirst(Filt.4$m)[, 3],col=2)みたいに、推定された状態の値の１列目と３列目を足し合わせた値を表示さ せています。
　じっさいのところ、推定されたFilt.4$mには１３列のデータがあります。そ のうち、「最終結果」にはどの列を使えばよいのかということは、モデルのFを見ればわかります。

> DLM.4$F
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
[1,]    1    0    1    0    0    0    0    0    0     0     0     0     0

これで、１列目と３列目だけが１になっているのがわかると思います。このFが何なのかを説明するには行列の掛け算の式を出さないといけない（実はかなり簡 単です。でも、このサイトでは行列を載せにくい。。。）ので省略しますが、とりあえずこれを見れば何を使えばよいかがわかります。
　ちなみに、１列目はローカルレベル＋トレンドの値、２列目はトレンドの値、３列目はその時期における季節成分です。だから１，３列目だけで十分だった と。
　４列目以降はその時期以外の季節成分たちです。補欠みたいなものですね（？）

　季節とかトレンドとかを別々のモデルで推定したので、結果も別々に表示させることができます。たとえばスムージングされた結果を使えば。

par(mfrow=c(3, 1))
#　元データ
plot(data, col=1, type="o", main="data")
#　レベル＋トレンド成分
plot(dropFirst(Smooth.4$s)[, 1], col=4, main="level+trend")
#季節成分
plot(dropFirst(Smooth.4$s)[, 3], col=4, main="seasonal")
par(mfrow=c(1, 1))

これを見れば、右肩上がりのトレンドがあること、そして、季節的な変動パタンが一目でわかります。

ｄｌｍな予測

次は予測をやってみます。１９５９年以降を切り落として、２年先までを予測することにします。

今度はシステムノイズも全部推定したモデルを作ります。
「全部」とは、ローカルレベルで説明した状態の変動と、トレンドの変動と、季節成分の変動です。
　今回は推定すべきパラメタが多いので、計算には多少の時間がかかります。また、最適化には多段階最適化法を使っていますが、それでも初期値の感度が大き くて、値を変えると別の結果が出たり、永遠に計算が終わらなかったりすることもあります。あくまでも暫定的な結果ということで見てください。

#　1959年以降を切り落とす
test.data <- window(data,end=c(1958,12))

#　Step1
#　モデル作成のための関数を作る
build.5<-function(theta){
    dlmModPoly(order=2,dV=exp(theta[1]),dW=c(exp(theta[2]),exp(theta[3])) )+
    dlmModSeas(fr=12,dW=c(theta[4],rep(0,10)),dV=0)
    }

#　Step2
#　MLEでパラメタ推定。
fit.5<-dlmMLE(test.data,parm=dlmMLE(test.data,parm=c(0,1,1,1),build.5,method="Nelder-Mead")$par,
　　　　　build.5,method="BFGS")

# 推定されたパラメタを使ってモデルを作り直す
DLM.5<-build.5(fit.5$par)

# Step3
# カルマンフィルター
Filt.5<-dlmFilter(test.data,DLM.5)

# Step4
# スムージング
Smooth.5<-dlmSmooth(Filt.5)

　なぜかフィルタリングの結果とスムージングの結果がほとんど一緒（わずかに違う値でした）でしたが、とりえあず気にしないで先に進んでいきます。

今回も、前に作った確定的なモデルと同じようにトレンドと季節に分けて表示させてやります。

par(mfrow=c(3,1))
#　元データ
plot(test.data,col=1,type="o", main="data")
#　レベル＋トレンド成分
plot(dropFirst(Smooth.5$s)[,1],col=4, main="level+trend")
#　季節成分
plot(dropFirst(Smooth.5$s)[,3],col=4, main="seasonal")
par(mfrow=c(1,1))

トレンドも季節もちょっとぶれてることがわかると思います。真ん中のトレンドのグラフは、季節調整済みデータとして使うこともできるでしょう。

で、次は本命の予測です。
dlmForecast()関数を使えば割と簡単に計算できます。

#　予測
Fore <- dlmForecast(Filt.5,nAhead=24,sampleNew=5)

#　予測の答え合わせ
plot(data, type="o")
lines(dropFirst(Smooth.5$s)[,1]+dropFirst(Smooth.5$s)[,3],col=4)
lines(Fore$f,col=2,lwd=2)
legend("bottomright",pch=c(1,NA),col=c(1,2),lwd=c(1,2),legend=c("実測値","予 測値"))

　ちょっと過小評価してますが、まぁまぁあっていそうです。
　ここで、dlmForecast(Filt.5,nAhead=24,sampleNew=5)と いうように書きましたが、nAhead=24で予測する期間を設定します。あとのsampleNew=5は、 予測値から乱数を使って出された予測値です。何のことかというと、予測値はあくまでも「予測された状態」です。なので、これに観測誤差なんかが加わって観 測値が出てくることになります。ということは、実際に観測値が出るとしたらどんな値が出るのか気になりますよね。それをシミュレーションしたんです。
　表示させるにはたとえばこんな風にすればよいです。

Line <- function(x){
    lines(x, col=8, type="o")
    }

plot(window(test.data,start=c(1957,1)), xlim=c(1957,1961),ylim=c(5.7,6.5), type="o")
lapply(Fore$newObs, Line)
lines(window(data,start=c(1959,1)), col=1,lwd=2)
lines(Fore$f,col=2)
legend("topleft",pch=c(NA,NA,1),col=c(1,2,8),lwd=c(2,1,1),legend=c("実測値 ","予測値","乱数を使った予測値"))

これをみると、予測結果からあり得なさそうな実測値が得られたわけではなさそうだということがわかります。これを使うと確率予報なんかも簡単にでき ますね。シミュレーション結果を１００個作って、そのうちの２０個が１００を超えたら、「１００を超える確率は２０％です」と予測してやればよいわけで す。近似的な確率ですが、予測の幅が広がります。



　バグや誤り等ございましたら、掲示板などにご一報ください。

参考文献

[1] J.J.Fコマンダー、S.Jクープマン　和合肇 訳：　状態空間時系列分析入門，２００８

[2] G,Petris, S Petrone, P Campagnoli: Dynamic Linear Models with R, ２００９

前のページ　⇒　ローカルレベルモデル　へ