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

前回やっていたちょっとエンジニアぽい試行（http://hsuetsugu.hatenablog.com/entry/2014/08/21/170046）がうまくいかなくなったので少し気分を変えて、分析モデルっぽいことで確認しておきたかったことをやってみます。

動機

小売や卸の日次・週次販売量予測をする際に、dlmモデルは一般的な時系列モデル（ARIMAなど）と比べてモデルの解釈が可能であり、かつ欠損値などに強い、という話を聞いたことがあります。これまで業務でdlmを使ったことはないのですが、販売量予測をする上で例えば店舗欠品のような無視すべきデータ*1のせいで、欠品後の予測値がぶれてしまうなどという経験は多くしてきたので、ずっと気になっていました。
また、価格を説明変数としたARIMAモデルなんかも使ったりしてきましたが、いまいちきれいなモデルが組めておらず、その辺がdlmをつかうとうまくできそうだ、という話も聞いていたので、この機会に試してみます。

dlm(Dynamic Linear Model)の説明

私自身理解できていませんが、入り口として下記のサイトが最もわかりやすかったです。
Logics of Blue
モデルの説明などは上記サイトなどをみていただくほうがいいと思いますが、ここでは実際にdlmとARIMAでのモデルフィッティングを通し、店舗欠品などがあった場合に対する予測値の違いなどをR言語にて検証していきたいと思います。

準備

ライブラリインストール

ARIMAの予測にforecast、dlmにはその名の通りdlmを使用します。内部処理にsqldfも使うかもしれません。

install.packages('forecast')
install.packages('dlm')
install.packages('sqldf')

dlmパッケージの作法の確認

dlmModPoly(dlmModPoly {dlm} | inside-R | A Community Site for R)

Usage

dlmModPoly(order = 2, dV = 1, dW = c(rep(0, order - 1), 1),
           m0 = rep(0, order), C0 = 1e+07 * diag(nrow = order))

Arguments

order:order of the polynomial model. The default corresponds to a stochastic linear trend.
dV:variance of the observation noise.
dW:diagonal elements of the variance matrix of the system noise.
m0:m0, the expected value of the pre-sample state vector.
C0:C0, the variance matrix of the pre-sample state vector.

dlmModSeas(dlmModSeas {dlm} | inside-R | A Community Site for R)

Usage

dlmModSeas(frequency, dV = 1, dW = c(1, rep(0, frequency - 2)),
           m0 = rep(0, frequency - 1),
           C0 = 1e+07 * diag(nrow = frequency - 1))

Arguments

frequency:how many seasons?
dV:variance of the observation noise.
dW:diagonal elements of the variance matrix of the system noise.
m0:m0, the expected value of the pre-sample state vector.
C0:C0, the variance matrix of the pre-sample state vector.

データ準備

あまりいい日次販売データがないのですが、今回は、Rにもともとある"BJsales"というデータセットを用いることにします。こんなデータです。

BJsales
Time Series:
Start = 1 
End = 150 
Frequency = 1 
  [1] 200.1 199.5 199.4 198.9 199.0 200.2 198.6 200.0 200.3 201.2 201.6 201.5 201.5 203.5
 [15] 204.9 207.1 210.5 210.5 209.8 208.8 209.5 213.2 213.7 215.1 218.7 219.8 220.5 223.8
 [29] 222.8 223.8 221.7 222.3 220.8 219.4 220.1 220.6 218.9 217.8 217.7 215.0 215.3 215.9
 [43] 216.7 216.7 217.7 218.7 222.9 224.9 222.2 220.7 220.0 218.7 217.0 215.9 215.8 214.1
 [57] 212.3 213.9 214.6 213.6 212.1 211.4 213.1 212.9 213.3 211.5 212.3 213.0 211.0 210.7
 [71] 210.1 211.4 210.0 209.7 208.8 208.8 208.8 210.6 211.9 212.8 212.5 214.8 215.3 217.5
 [85] 218.8 220.7 222.2 226.7 228.4 233.2 235.7 237.1 240.6 243.8 245.3 246.0 246.3 247.7
 [99] 247.6 247.8 249.4 249.0 249.9 250.5 251.5 249.0 247.6 248.8 250.4 250.7 253.0 253.7
[113] 255.0 256.2 256.0 257.4 260.4 260.0 261.3 260.4 261.6 260.8 259.8 259.0 258.9 257.4
[127] 257.7 257.9 257.4 257.3 257.6 258.9 257.8 257.7 257.2 257.5 256.8 257.5 257.0 257.6
[141] 257.3 257.5 259.6 261.1 262.9 263.3 262.8 261.8 262.2 262.7

予測の前提

直近98日間(14週間)のデータを用いて翌週7日間の日次販売量を予測していくことにします。今回用いるデータは150日間なので、後ろ7週間を予測対象期間として毎週ローリングしていくこととします。

予測結果

パラメータの推定はauto.arimaに完全に任せて、上記の前提で7週間分の予測をまわしました。結果は下記のような感じです。結構なトレンドのなかで、ARIMAもdlmもそんなに悪くはない感じがします。
f:id:hsuetsugu:20140822163144p:plain

コードは下記の通りです。ちなみにplotするときのlegendのlocator(1)というのを初めて知りました・・！便利です。

data <- BJsales
# plot(data, type='l', lwd=2)

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

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

for(i in 1:targetduration){
  data.this <- data[(7*(i-1)+1):(7*(prevduration + i - 1))]

  ### 1. ARIMAモデル最適パラメーター探索
  fit <- try(auto.arima(data.this, ic="aic", trace=F, stepwise=F, start.p=0, start.q=0, start.P=0, start.Q=0), silent=T)
  if(as.vector(attributes(fit)[1])=="try-error") next
  pred.arima <- try(forecast(fit, level = c(95), h = 7), silent=T)
  if(as.vector(attributes(pred.arima)[1])=="try-error") next    
  x.arima <- append(x.arima,list(pred.arima))  
  
  ### 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)
  }
  fit0 <- try(dlmMLE(data.this, parm=c(0,1,1), bmodel), TRUE)
  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))   
  
}

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

#予測値の格納(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)')
lines(frcst.arima, col=2, lwd=2)
lines(frcst.dlm, col=4, lwd=2)
legend(locator(1), c('actual','arima','dlm'), pch=c(1,1,1), col=c(1,2,4),cex=1)
grid()