最終更新:2017年6月1日

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

●状態空間モデル関連のページ
なぜ状態空間モデルを使うのか
状態空間モデル:状態空間モデルのことはじめ
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になったとします。
すると、2001年の状態は

10(元の値) + 3(状態の変動) + 4-2(トレンド) = 15   となります。

ノイズの大きさは常に±7なので、2001年の観測値は、8~22の範囲内に収まるだろうと考えられます。

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

10(2000年の状態の値) + 0 + 4(2000年のトレンドの値) = 14

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

 

4.季節成分

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

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

 

5.dlmで季節とトレンド

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

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

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

データのプロット

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

dlmパッケージを落としてから

library(dlm)

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

# レベルも傾きも、季節変動も、システムノイズ0 観測誤差のみ
# 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ですが、これが1ならば単なるローカルレベルモデルで、2にするとトレンド付きになります。ちなみにデフォルトは2なので、いちいち書かなくても実は問題ありません。
dVは観測誤差で、この大きさは推定しています。でも、状態そのものは確率的に変動しないというモデルを今回は組んでいるので、dWは0を指定しています。ここで、dWに二つも0を入れていますが、一つ目はローカルレベルモデルで説明した状態の変動の大きさを表し、二つ目はトレンドの変動の大きさを指定しています。ですので、片方だけ0とか言ったモデルももちろん推定できます。

紫の部分は季節を表したところですが、これにはdlmModSeas()を使います。fr=12とは、一年の間に12こデータがありますよという 指示。たとえば春夏秋冬4つの区分しかなければfr=4とします。
dWにやけにたくさん0があります(最初の0に追加してrep()でさらに0を10個生産している)が、これは最後の10個の0はほぼ確定です。気にせ ず0を入れればよいです。
dWの一番最初の値はまた別です。季節が確率的変動をする場合は、この大きさも推定することになります。逆に言えば、季節が確率的変動をしない場合でも、最後の0を10個くっつけるところは変化しないことに気を付けてください。
観測誤差のほうはdlmModPoly()で推定されているので、ここには0を突っ込んでおきます。

推定された観測誤差は

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

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

グラフを書きます

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

legend("bottomright", pch=c(1,NA,NA),
col=c(1,2,4), lwd=c(1,2,2), legend=c("data","Filter","Smooth"))

プロセスエラー無し

lines(dropFirst(Filt.4$m)[, 1] + dropFirst(Filt.4$m)[, 3], col=2, lwd=2)のように、推定された状態の値の1列目と3列目を足し合わせた値を表示させています。

じっさいのところ、推定されたFilt.4$mには13列のデータがあります。そ のうち、「最終結果」にはどの列を使えばよいのかということは、モデルのFを見ればわかります。

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

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

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

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))

プロセスエラー無しその二

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

 

6.dlmな予測

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

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

# 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("実測値","予測値","乱数を使った予測値"))

予測その二

 

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

バグや誤り等ございましたら、メールなどでご一報いただけると幸いです。

 

参考文献


Dynamic Linear Models with R (Use R!)

 
私がdlm関連の記事を書くときは、いつもこの本を参考にしています。
dlmパッケージを使いこなしたい方はぜひどうぞ。
 


Rによるベイジアン動的線形モデル (統計ライブラリー)

 
Dynamic Linear Models with Rの日本語訳です。
 

状態空間時系列分析入門

 
このサイトでいつも紹介している、状態空間モデルの入門書です。
状態空間モデルについて知りたければ、まずはこの本から始めるのが良いと思います。最も読みやすい状態空間モデルの入門書です。
 

状態空間モデル関連のほかの記事はこちらから見ることができます



スポンサードリンク