(35連休)10日目：ARIMAモデルとdlm（動的線形モデル）による販売量予測モデルの検証③

前回（http://hsuetsugu.hatenablog.com/entry/2014/08/23/125238）の続きです。

前回のおさらいと訂正

前回、RのBJsalesというデータに対して意図的に欠損値をいれて、dlmとARIMAでの予測結果を比較しました。動機は、繰り返しになりますが、「状態」を推定するdlmでは欠損値に強い、というのを検証することです。
前回計算した結果は下記ですが、設定したパラメータに問題があることがわかりましたので今回はその訂正の記事を書きます。

訂正①：dlmModPolyの部分

推定しているローカルレベルモデルのdV（観測誤差の大きさ）とdW（状態の変動の大きさ）について、初期値を与えているのですが、ここの初期値の影響を受けていて、正しい最適化がでていない、ということになっていました。dlmMLEの結果を受けているfit0について、fit0$mにてカルマンフィルターで推定された状態をみてみると、原系列と同じになっていてようやく気づいたのですが・・・。
そこで、下記のように訂正しました。

  fit0 <- dlmMLE(
    data.this,
    parm=dlmMLE(data.this, parm=c(0,0), bmodel,method='Nelder-Mead')$par,bmodel, 
    method='BFGS'
  )

これは、「多段階最適化」というものらしく、

dlmMLEを使って普通にパラメタを求めた後、その推定されたパラメタを初期値にしてもう一度最適化してやる方法

ということです。（ローカルレベルモデル | 状態空間モデル | Logics of Blue）

訂正②：dlmModSeasの部分

今回のデータBJsalesではあまりそういう傾向はなさそうではあるのですが、だいたいの販売データは強い曜日変動があります。なので、dlmModSeasのfr=7で曜日変動を表したかったのですが、これを指定するとどうしても欠損値に予測値がふれてしまうので、一旦季節変動の項自体なくしました（モデルは上述のコードの通り）。

結果

ARIMAモデルの結果は省いて、各週で推定したカルマンフィルターで推定された状態と、翌７日間の予測値をプロットしたものが下記です。直前のカルマンフィルタで推定された状態の値そのものになってそうですが、そういうことなのでしょう。どうであれ、結構な欠損値があるにも関わらずある程度安定した予測値を算出してくれたのでまずはよかったです。引き続き、正しく理解できるよう努めます。

変更点は上述の通りですが、一応コードをはっておきます。

data <- BJsales
# data[50:51]<-0
  
# plot(data, type='l', lwd=2)

targetduration <- 7 # 7w予測する
prevduration <- 14 #直近14wを使う

### 初期化
x.dlm <- as.list(NULL)
x.dlm.m <- as.list(NULL)

for(i in 1:targetduration){
  # i <- 1
  data.this <- data[(7*(i-1)+1):(7*(prevduration + i - 1))]
  
  ### 2. 状態空間モデル
  #### トレンドは推定しない場合 #### 
  bmodel <- function (theta) {
    # dlmModPoly (order=1, dV=exp(theta[1]), dW=exp(theta[2])) + dlmModSeas(fr=7, dW=c(theta[3], rep(0,5)), dV=0)
    # dlmModPoly (order=1, dV=exp(theta[1]), dW=exp(theta[2])) + dlmModSeas(fr=98, dW=c(theta[3], rep(0,96)), dV=0)
    dlmModPoly (order=1, dV=exp(theta[1]), dW=exp(theta[2])) # ローカルレベルモデル
  }
  fit0 <- dlmMLE(
    data.this,
    # parm=dlmMLE(data.this, parm=c(1,1,1), bmodel,method='Nelder-Mead')$par,bmodel,
    parm=dlmMLE(data.this, parm=c(0,0), bmodel,method='Nelder-Mead')$par,bmodel, #ローカルレベルモデル
    method='BFGS'
  )
  if(is.null(attributes(fit0)$class)){
    ssmodel <- bmodel(fit0$par)
    fit <- dlmFilter(data.this, ssmodel)
    pred.dlm <- dlmForecast(fit, nAhead=7)
  } else {
    pred.dlm <- list(f=rep(NA, 7))
  }
  x.dlm <- append(x.dlm,list(pred.dlm))   
  x.dlm.m <- append(x.dlm.m,list(fit$m))
}

#予測値の格納(dlm)
frcst.dlm <- as.vector(NULL)
for(i in 1:targetduration){
  frcst.dlm <- append(frcst.dlm,x.dlm[[i]]$f)
}
frcst.dlm <- append(rep(NA,prevduration*7),frcst.dlm)

# 結果描画
plot(data, ylim=c(200,270),type='l', lwd=2,main='result (arima / dlm)')
for (j in 1: targetduration){ #推定された状態の表示
  lines(c(0:(7*(j-1)),dropFirst(x.dlm.m[[j]])), col=j+1, lwd=2)  
}
lines(frcst.dlm, col=4, lwd=2)
legend(locator(1), c('actual','dlm'), pch=c(1,1), col=c(1,4),cex=1)
grid()