状態空間モデルのマーケティングへの応用（サンプルコード）

#######参考プログラム#######

#データの読み込み:test.csvのファイルレイアウト

#ID_W:週ID

#end:エンド陳列の有無

#logPI:点数PIの対数

#AP:売価

#RP095:参照価格(繰越パラメータ0.95：RP_t=0.95*RP_{t-1}+(1-0.95)*AP_{t-1})

#

#（注）下記プログラムでは，簡便のため定常AR成分p_tは除いている

#（注）本文中のデータとは異なっている

#

library(dlm)

Dataset <- read.table("test.csv", header=TRUE,sep=",", na.strings="NA", dec=".", strip.white=TRUE)

#フルモデル（RP095）

n<-nrow(Dataset)

y<-Dataset$logPI

E<-Dataset$end

#価格変数の加工

#実売価と参照価格の差の算定

#

z_2<-0

z_3<-0

v_wk<-Dataset$AP-Dataset$RP095

#価格変数の算定（ロス変数）

for (i in 1:n){

if (v_wk[i]>=0) {

z_2[i] <- (Dataset$AP[i]-Dataset$RP095[i])/100

} else {

z_2[i] <- 0

}

#価格変数の算定（ゲイン変数）

if (v_wk[i] < 0) {

z_3[i] <- (Dataset$RP095[i]-Dataset$AP[i])/100

} else {

z_3[i] <- 0

}

}

#

#価格変数の交互作用変数の設定

z_2E<-z_2*E

z_3E<-z_3*E

#

#デザイン行列の設定

x<-cbind(z_2,z_2E,z_3,z_3E)

#

#モデルの設定

buildModel <- function(params){

#dlmModRegには説明変数を入れる

mod1 <- dlmModReg(x, dV = exp(params[1]), dW = exp(params[2:6]))

return(mod1)

}

#最尤法

outMLE1 <- dlmMLE(y,rep(0, 6),buildModel,method = "L-BFGS-B",hessian=T,control = list(maxit = 1000000,trace=1))

#最尤法の結果

outMLE1$convergence

avar<-solve(outMLE1$hessian)

outMLE1$par

outMLE1$value

AIC<-2*outMLE1$value+2*6

AIC

exp(outMLE1$par)

sqrt(diag(avar))

#最尤推定結果に基づくフィルタリング

#

mod1 <- buildModel(outMLE1$par)

mod.filt1<-dlmFilter(y,mod1)

#固定区間平滑化

mod.smooth1 <- dlmSmooth(mod.filt1)

a1<-dropFirst(mod.smooth1$s[,1])

a2<-dropFirst(mod.smooth1$s[,2])

a3<-dropFirst(mod.smooth1$s[,3])

a4<-dropFirst(mod.smooth1$s[,4])

a5<-dropFirst(mod.smooth1$s[,5])

#平滑化のグラフ

par(mfrow=c(3,2))

plot(a1,type="l",xlab="weeks",ylab="beta_0" )

plot(a2,type="l",xlab="weeks",ylab="beta_1")

plot(a3,type="l",xlab="weeks",ylab="beta_2")

plot(a4,type="l",xlab="weeks",ylab="beta_3")

plot(a5,type="l",xlab="weeks",ylab="beta_4")