#######参考プログラム#######
#データの読み込み: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")