最終更新:2017年6月1日dlmパッケージを使って、ちょっと複雑な正規線形状態空間モデル(動的線形モデル)を推定します。今回は、トレンド成分と季節成分の入った状態空間モデルを推定し、予測を行うところまでを行います。:Rで正規線形状態空間モデルを当てはめるローカルレベルモデル:dlmパッケージを使ってローカルレベルモデルを当てはめる季節とトレンド:dlmパッケージを使って季節成分とトレンドの入ったモデルを作るdlmによる時変形数モデル:dlmによる「時間によって係数が変化する回帰モデル」の作成Pythonによる状態空間モデル:R言語ではなくPythonを使いたい方はこちらをどうぞスポンサードリンク目次1.ローカルレベルモデルを発展させる2.トレンド3.ローカル線形トレンドモデル4.季節成分5.dlmで季節とトレンド6.dlmな予測1.ローカルレベルモデルを発展させるこちらで書いたように、一番単純な動的線形モデルであるローカルレベルモデルはデータから正しそうな状態を推定するということはある程度可能でしたが、予測には何の役にも立ちませんでした。これは、モデルを作るのに使った「データに対する知識」があまりにも少なかったことが原因です。「このデータは上昇(下降)トレンドがある」とか、「毎年8月には高い値を示す」といった季節変動といったの知識を取り込むとより良い予測が可能です。焦点を当てて、これらの要因を動的線形モデルに付け加える方法について記します。2.トレンドトレンドは大きく分けて2つあります。一つは、確定的トレンドもう一つが確率的トレンドです。確定的トレンドは、回帰分析を思い出してもらえるとよいです。散布図に直線を引っ張るアレです。横軸に時間を取って、縦軸にたとえばビールの売り上げを置いたら、それこそ右肩上がりに一定のスピードで増加していく、というのが確定的トレンド。図で示すとこんな感じです。上のグラフの場合は、時間が1増えるごとに、ビールの売り上げは2ずつ増えていっていることがわかると思います。一定して、時間が1増えると、売り上げが2増える、というのが、確定的トレンドです。確率的トレンドは、その名の通り、確率的に増加(減少)分が変化します。この変化のパターンもランダムウォークです。一期前のトレンドと似たようなトレンドになるけれども、ちょっと値が変わるということです。3.ローカル線形トレンドモデルトレンドとローカルレベルモデルを合体させたもののイメージをつかんでいただくために、例を使って説明します。2000年におけるビールの売り上げは10本。その時のトレンドは、+4でした。順当にいけば2001年は14本売れることになります。最尤法をした結果、観測誤差は±7の範囲を持つことがわかり、状態の変動の大きさは±4であることが推定できたとします。トレンドの変動の大きさは±2でした。2000年の状態は10でした。2001年は、10から±4の範囲内で動く可能性があります。トレンドは4±2なので、2~6の範囲で変動します。サイコロを振ったら状態の変化は+3になったとします。トレンドの変化は-2になったとします。ノイズの大きさは常に±7なので、2001年の観測値は、8~22の範囲内に収まるだろうと考えられます。予測を行う際、サイコロは正負どちらに転ぶかわかりません。いう時は「変化しない」を前提として予測するんでした。14という風に予測されます。2001年のデータが手に入ったら、それをもとにして予測された状態を補正する、という流れは一緒です。今度は季節成分について説明します。季節成分の入れ方には三角関数(サインとかコサインとか)を使ったやり方と、ダミー変数を使ったやり方の2種類がありますが、ここでは後者のダミー変数を使う方法だけを説明します(dlmパッケージを使って三角関数な季節成分を推定することも可能です)。ダミー変数を使ったやり方とは、単純にひと月ごとに売り上げの大きさの補正をすることです。1月にはビールは売れないけど、8月にはよく売れる、といったような感じ。この季節成分にもシステムノイズを入れ込むことは可能です。ノイズがあるのが確率的季節成分、ないのが確定的季節です。1月には売れ行きが悪くて-5の補正、2月には若干ましになって-3の補正を受けたとします。2000年1月の売り上げは10でした。季節成分だけを見ると、2月は1月よりプラス2になっているので、2000年2月の売り上げは12となるはずです。これがシステムノイズがない場合。システムノイズがあるときは、ローカルレベルモデルなんかと同じく、季節ごとにさいころを投げてやります。2000年2月の値は、11~13ということになります。この季節のブレだけでなく、ローカルレベルモデルで表された状態の変動やトレンドの変動なんかも加えていくことが可能です。5.dlmで季節とトレンド今回使ったソースはこちらにおいてあります。今回は、季節成分とトレンド両方がありそうなデータとして、AirPassengersというデータを使います。データを対数変換したものを解析の対象として使います。lwd=2)まずは、状態の変化に一切確率的な変動がない、プロセスノイズが全部0なモデルを推定します。ローカルレベルモデルで示したよう、システムノイズが0な状態空間モデルは、普通の線形回帰をした場合と一緒の結果になります。このモデルは重回帰分析を使うことによっても推定可能です。暇な方は試してみてください。dlmパッケージを落としてからlibrary(dlm)あとは、以前に示したように4ステップで推定完成です。MLEでパラメタ推定。dlmSmooth(Filt.4)モデルの型を決める作業が一番大きく変わったと思います。dV=0)}において、緑の部分がローカル線形トレンドモデル(ローカルレベルモデル+トレンド)で、紫が季節の成分です。つなげば簡単に合体できます。緑の部分から説明します。太線になっているorder=2ですが、これが1ならば単なるローカルレベルモデルで、2にするとトレンド付きになります。ちなみにデフォルトは2なので、いちいち書かなくても実は問題ありません。dVは観測誤差で、この大きさは推定しています。状態そのものは確率的に変動しないというモデルを今回は組んでいるので、dWは0を指定しています。ここで、dWに二つも0を入れていますが、一つ目はローカルレベルモデルで説明した状態の変動の大きさを表し、二つ目はトレンドの変動の大きさを指定しています。片方だけ0とか言ったモデルももちろん推定できます。紫の部分は季節を表したところですが、これにはdlmModSeas()を使います。春夏秋冬4つの区分しかなければfr=4とします。dWにやけにたくさん0があります(最初の0に追加してrep()でさらに0を10個生産している)が、これは最後の10個の0はほぼ確定です。ず0を入れればよいです。dWの一番最初の値はまた別です。季節が確率的変動をする場合は、この大きさも推定することになります。逆に言えば、季節が確率的変動をしない場合でも、最後の0を10個くっつけるところは変化しないことに気を付けてください。観測誤差のほうはdlmModPoly()で推定されているので、ここには0を突っ込んでおきます。0.0035172となりました。ふつうにDLM.4と推定されたモデルの型を表示させるとすごくたくさんの数値が出てきますが、基本的な見方はローカルレベルモデルと一緒なので省略します。季節があるから、その分表示する量が多くなったんですね。lwd=2)のように、推定された状態の値の1列目と3列目を足し合わせた値を表示させています。じっさいのところ、推定されたFilt.4$mには13列のデータがあります。うち、「最終結果」にはどの列を使えばよいのかということは、モデルのFを見ればわかります。0これで、1列目と3列目だけが1になっているのがわかると思います。このFが何なのかを説明するには行列の掛け算の式を出さないといけない(実はかなり簡単です。このサイトでは行列を載せにくい)ので省略しますが、とりあえずこれを見れば何を使えばよいかがわかります。ちなみに、1列目はローカルレベル+トレンドの値、2列目はトレンドの値、3列目はその時期における季節成分です。4列目以降はその時期以外の季節成分たちです。補欠みたいなものですね(?)季節とかトレンドとかを別々のモデルで推定したので、結果も別々に表示させることができます。スムージングされた結果を使えば。1))これを見れば、右肩上がりのトレンドがあること、そして、季節的な変動パタンが一目でわかります。次は予測をやってみます。1959年以降を切り落として、2年先までを予測することにします。今度はシステムノイズも全部推定したモデルを作ります。「全部」とは、ローカルレベルで説明した状態の変動と、トレンドの変動と、季節成分の変動です。今回は推定すべきパラメタが多いので、計算には多少の時間がかかります。最適化には多段階最適化法を使っていますが、それでも初期値の感度が大きくて、値を変えると別の結果が出たり、永遠に計算が終わらなかったりすることもあります。暫定的な結果ということで見てください。MLEでパラメタ推定。dlmSmooth(Filt.5)なぜかフィルタリングの結果とスムージングの結果がほとんど一緒(わずかに違う値でした)でしたが、とりえあず気にしないで先に進んでいきます。今回も、前に作った確定的なモデルと同じようにトレンドと季節に分けて表示させてやります。main=”seasonal”)par(mfrow=c(1,1))トレンドも季節もちょっとぶれてることがわかると思います。真ん中のトレンドのグラフは、季節調整済みデータとして使うこともできるでしょう。dlmForecast()関数を使えば割と簡単に計算できます。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で予測する期間を設定します。予測値から乱数を使って出された予測値です。何のことかというと、予測値はあくまでも「予測された状態」です。これに観測誤差なんかが加わって観測値が出てくることになります。いうことは、実際に観測値が出るとしたらどんな値が出るのか気になりますよね。それをシミュレーションしたんです。legend=c(“実測値”,”予測値”,”乱数を使った予測値”))これをみると、予測結果からあり得なさそうな実測値が得られたわけではなさそうだということがわかります。