[R] bsts posterior distributions

[Greg Dropkin]

hello

I couldn't find previous posts to answer this, but please point me to any.

I am trying to understand bsts, starting with no regressors

Here is code which appears to mimic bsts, producing graphs similar to the
model plot, but gives a rather different posterior distribution for the
parameters. I may misunderstand SdPrior, or perhaps the differences are
with mcmc.

Thanks for any help

Greg

###

library(bsts)
library(Boom)
library(mcmc)

#simulate some data

y<-rep(NA,50)

y[1]=1
y[2]=1
s=2
set.seed(5)
for (k in 1:48)
{
y[2+k]=y[1+k]+0.1*y[k]+s*rnorm(1)
}
plot(1:50,y[1:50],main=paste("seed =",j))

#bsts model

ss<-AddLocalLevel(list(),y)
mod1<-bsts(y,state.specification=ss,niter=1000)
plot(mod1)

#reproduce the plot using mod1$state.contributions

par(mfrow=c(1,2))
plot(1:50,y,col="blue",ylim=c(-12,8),main="quantiles of
mod1$state.contributions")
for (i in 1:99)
{
qi<-qin<-rep(NA,50)
tj<-1:50
for (j in 1:50)
{
qi[j]<-quantile(mod1$state.con[501:1000,1,j],(i-0.5)/100)
qin[j]<-quantile(mod1$state.con[501:1000,1,j],(i+1-0.5)/100)
}
polygon(c(tj,rev(tj)),c(qi[1:50],rev(qin[1:50])),col=rgb(0,0,0,45*dnorm(i,mean=50,sd=18)),border=FALSE)
}
plot(mod1,ylim=c(-12,8),main="plot(mod1)")
par(mfrow=c(1,1))

#the state specification is based on sd(y)

sd(y)
ss

#the likelihood, using the kalman filter, as a function of the error sd's

kal<-function(par)
{
a<-par[1]
b<-par[2]

H=matrix(1,1,1)
F=matrix(1,1,1)
#1-dimensional state
N=50
dim(y)=c(1,N)

xe<-ye<-ze<-matrix(NA,1,N)
w<-rnorm(1)
xe[,1]<-1
ye[,1]<-H%*%xe[,1]
ze[,1]<-y[,1]-ye[,1]

P<-K<-array(data=NA,dim=c(N,1,1))
#P[1,,] initial guess
P[1,,]<-1
K[1,,]<-1
xe[1,1]<-1

for (i in 1:(N-1))
{
P[i+1,,]<-F%*%P[i,,]%*%t(F)+a^2
K[i+1,,]<-P[i+1,,]%*%t(H)%*%solve(b^2+H%*%P[i+1,,]%*%t(H))
xe[1,i+1]<-F%*%xe[,i]+K[i+1,,]%*%(y[,i+1]-H%*%F%*%xe[,i])
P[i+1,,]<-(diag(1,1)-K[i+1,,]%*%H)%*%P[i+1,,]
}

-1/2*log(b^2)-1/2*sum(((y[1,]-xe[1,])/b)^2)
}

#independent priors on a and b
#bsts uses sd(y) to set both priors.
#for a (ss[[1]]$sigma.prior) it uses SdPrior with
#$prior.guess 0.04600655, $prior.df 0.01, $initial.value 0.04600655,
$upper.limit 4.600655
#for b (mod1$prior) it uses SdPrior with
#$prior.guess 4.600655, $prior.df 0.01, $initial.value 4.600655,
$upper.limit 5.520786

#try an inverse gamma prior (on a^2, then converted to a prior on a)

lpr1<-function(a)
{
v=0.01
ifelse((a<=0)|a>sd(y),-Inf,-(v/2+1)*log(a^2)-v*(sd(y)/100)^2/(2*a^2)+1/2*log(a^2))
}

lpr2<-function(b)
{
v=0.01
ifelse((b<=0)|b>1.2*sd(y),-Inf,-(v/2+1)*log(b^2)-v*(sd(y))^2/(2*b^2)+1/2*log(b^2))
}

lpost1<-function(par)
{
a<-par[1]
b<-par[2]
lpr1(a)+lpr2(b)+kal(par)
}

nlpost1<-function(par) -lpost1(par)

pop<-optim(c(2,1),nlpost1)
pop<-optim(pop$par,nlpost1,hessian=T)
pop

lpost1m<-function(par) lpost1(par+pop$par)
nlpost1m<-function(par) -lpost1m(par)
pop1m<-optim(c(0,0),nlpost1m)
pop1m<-optim(pop1m$par,nlpost1m,hessian=T)
pop1m

sc<-chol(pop1m$hess)
sd<-diag(sc)
nb=5000
#with batches, spacing, ~5 min run
outm<-metrop(lpost1m,1e-2*sd,nb,blen=5,nspac=5)
samm<-outm$batch
dim(samm)
acf(samm)
samm<-thin(samm,5)
dim(samm)
acf(samm)
par(mfrow=c(2,1))
plot(samm[501:1000,1],type="l")
plot(samm[501:1000,2],type="l")
par(mfrow=c(1,1))

#plot kalman using samm from lpost1m

ax<-samm[501:1000,1]+pop$par[1]
bx<-samm[501:1000,2]+pop$par[2]

state<-matrix(NA,500,50)

for (j in 1:500)
{
a<-ax[j]
b<-bx[j]

H=matrix(1,1,1)
F=matrix(1,1,1)
N=50
dim(y)=c(1,N)

xe<-ye<-ze<-matrix(NA,1,N)
xe[,1]<-1
ye[,1]<-H%*%xe[,1]
ze[,1]<-y[,1]-ye[,1]

P<-K<-array(data=NA,dim=c(N,1,1))
P[1,,]<-0
for (i in 1:(N-1))
{
P[i+1,,]<-F%*%P[i,,]%*%t(F)+a^2
K[i+1,,]<-P[i+1,,]%*%t(H)%*%solve(b^2+H%*%P[i+1,,]%*%t(H))
P[i+1,,]<-(diag(1,1)-K[i+1,,]%*%H)%*%P[i+1,,]
xe[1,i+1]<-F%*%xe[,i]+K[i+1,,]%*%(y[,i+1]-H%*%F%*%xe[,i])
}

P[1,,]<-P[2,,]
state[j,]<-rnorm(N,xe[1,1:N],P[1:N,,]^0.5)
}

par(mfrow=c(1,2))
plot(1:N,y,col="blue",ylim=c(-12,8))
for (i in 1:99)
{
qi<-qin<-rep(NA,N)
tj<-1:N
for (k in 1:N)
{
qi[k]<-quantile(state[,k],(i-0.5)/100)
qin[k]<-quantile(state[,k],(i+1-0.5)/100)
}
polygon(c(tj,rev(tj)),c(qi[1:N],rev(qin[1:N])),col=rgb(0,0,0,45*dnorm(i,mean=50,sd=18)),border=FALSE)
}
plot(mod1,ylim=c(-12,8))
par(mfrow=c(1,1))

#very plausible, but

mean(ax)
mean(mod1$sigma.level[501:1000])

mean(bx)
mean(mod1$sigma.obs[501:1000])

par(mfrow=c(1,2))
qqplot(ax,mod1$sigma.level[501:1000])
lines(ax,ax,col="red")
qqplot(bx,mod1$sigma.obs[501:1000])
lines(bx,bx,col="red")
par(mfrow=c(1,1))

#ax not sampling sigma.level
#bx not quite sampling sigma.obs