First order DLM

Single-move vs block-move

1 Local level model set up

The following exercise is inspired in Gamerman, Reis and Salazar (2006) Comparison of sampling schemes for dynamic linear models. International Statistical Review, 74, 203-214.

The purpose of the exercise is to compare two MCMC schemes to sample states and fixed parameters in the first order dynamic linear model, also known as the local level model that we have already discussed in class. More specifically, the local level model assumes that the time series data, \(\mbox{data}=\{y_1,\ldots,y_n\}\), can be described as follows, for \(t=1,\ldots,n\) \[\begin{eqnarray*} y_t &=& \beta_t + \nu_t \qquad \nu_t \sim N(0,V)\\ \beta_t &=& \beta_{t-1} + \omega_t \qquad \omega_t \sim N(0,W), \end{eqnarray*}\] with \(\beta_0 \sim N(m_0,C_0)\) and where \(W/V\) is commonly known as the signal to noise ratio (SNR). Here \(m_0\) and \(C_0\) are known hyperparameters.

The model is completed with prior distributions for the error variances, \(V\) and \(W\). Here, we will assume that \[\begin{eqnarray*} V &\sim& IG(a_\nu,b_\nu)\\ W &\sim& IG(a_\omega,b_\omega), \end{eqnarray*}\] \(a_\nu,b_\nu,a_\omega,b_\omega\) known hyperparameters.

2 Simulating one dataset

For simplication, we created the data simulation function sample.data.

sample.data = function(n,W,V,beta0){
  sdW     = sqrt(W)
  sdV     = sqrt(V)
  beta    = rep(0,n)
  y       = rep(0,n)
  beta[1] = rnorm(1,beta0,sdW)
  y[1]    = rnorm(1,beta[1],sdV)
  for (t in 2:n){
    beta[t] = rnorm(1,beta[t-1],sdW)  
    y[t]    = rnorm(1,beta[t],sdV)
  }
  return(y)
}

We simulated \(n=100\) observations, starting with \(\beta_0=0\) and fixing \(V=1.0\) and \(W=0.25\) so the SNR is 0.25.

set.seed(17389)
n     = 100
V     = 1.0
W     = 0.25
beta0 = 0.0
y     = sample.data(n,W,V,beta0)
ts.plot(y)

3 Block-move

The following routine mcmc.joint samples \((\beta_1,\ldots,\beta_n)\), \(V\) and \(W\) iteratively from their full conditional posterior distributions, ie. \[\begin{eqnarray*} p(\beta_1,\ldots,\beta_n|V,W,\mbox{data})\\ p(V|\beta_1,\ldots,\beta_n,W,\mbox{data})\\ p(W|\beta_1,\ldots,\beta_n,V,\mbox{data}), \end{eqnarray*}\] and returns an \(M \times (n+2)\) matrix with \(M\) draws of the \((n+2)\)-dimensional vector \((\beta_1,\ldots,\beta_n,V,W)\).

The full conditional distributions of \(V\) and \(W\) are both Inverse-Gamma distributions, while the full conditional distribution of \((\beta_1,\ldots,\beta_n)\) is \(n\)-variate normal distribution and its joint sampling is facilitated by the kalman recursions, both filtering and smoothing.

Below, \(a_1=m_0\) and \(R_1=C_0+W\), while \(M_0\) and \(M\) are MCMC draws, with the first \(M_0\) draws discarded (burn-in or warm-up) and the following \(M\) draws stores for posterior inference. Initial values for \(V\) and \(W\) are the two last arguments of the function mcmc.joint.

mcmc.joint = function(y,a1,R1,av,bv,aw,bw,M0,M,V,W){
  n      = length(y)
  V.draw = V
  W.draw = W
  niter  = M0+M
  draws  = matrix(0,niter,n+2)
  for (iter in 1:niter){
    beta.draw = ffbs(y,rep(1,n),0.0,V.draw,1.0,0.0,W.draw,a1,R1,nd=1)
    par2v     = bv+sum((y-beta.draw)^2)/2
    par2w     = bw+sum((beta.draw[2:n]-beta.draw[1:(n-1)])^2)/2
    V.draw    = 1/rgamma(1,av+n/2,par2v)
    W.draw    = 1/rgamma(1,aw+(n-1)/2,par2w)
    draws[iter,] = c(beta.draw,V.draw,W.draw)
  }
  return(draws[(M0+1):niter,])
}

3.1 Sampling from \(p(\beta_1,\ldots,\beta_n|V,W,\mbox{data})\)

The following routine ffbs was actually written to be able to sample the states jointly via the forward filtering backward sampling (FFBS) scheme for a more general class of normal dynamic linear model (NDLM):

\[\begin{eqnarray*} y_t &=& F_t^T \theta_t + \nu_t \qquad \nu_t \sim N(0,V)\\ \theta_t &=& \gamma + G \theta_{t-1} + \omega_t \qquad \omega_t \sim N(0,W), \end{eqnarray*}\] where \(F_t\) is a vector of observed components, with \(\theta_t\) a state vector of the same dimensions. If \(F_t\) are the usual definition of regressors in linear regression models, then \(\theta_t\) are the time-varying regression coefficients. For simple local level model, we have that \(F_t=1, \ \ \forall t\), \(\gamma=0\), \(G=1\) and \(\theta_t=\beta_t\). Another special case, that we have studied in this class, is the AR(1) plus noise model: \[\begin{eqnarray*} y_t &=& \theta_t + \nu_t \qquad\qquad \nu_t \sim N(0,V)\\ \theta_t &=& \gamma + \phi \theta_{t-1} + \omega_t \qquad \omega_t \sim N(0,W). \end{eqnarray*}\]

We use the function ffbs to draw \(\beta_1,\ldots,\beta_n\) in the MCMC function mcmc.joint.

ffbs = function(y,Ft,alpha,V,G,gamma,W,a1,R1,nd=1){
  n = length(y)
  if (length(Ft)==1){Ft=rep(Ft,n)}
  if (length(alpha)==1){alpha=rep(alpha,n)}
  if (length(V)==1){V=rep(V,n)}
  a = rep(0,n)
  R = rep(0,n)
  m = rep(0,n)
  C = rep(0,n)
  B = rep(0,n-1)
  H = rep(0,n-1)
  # time t=1
  a[1] = a1
  R[1] = R1
  f    = alpha[1]+Ft[1]*a[1]
  Q    = R[1]*Ft[1]**2+V[1]
  A    = R[1]*Ft[1]/Q
  m[1] = a[1]+A*(y[1]-f)
  C[1] = R[1]-Q*A**2
  # forward filtering
  for (t in 2:n){
    a[t] = gamma + G*m[t-1]
    R[t] = C[t-1]*G**2 + W
    f    = alpha[t]+Ft[t]*a[t]
    Q    = R[t]*Ft[t]**2+V[t]
    A    = R[t]*Ft[t]/Q
    m[t] = a[t]+A*(y[t]-f)
    C[t] = R[t]-Q*A**2
    B[t-1] = C[t-1]*G/R[t]
    H[t-1] = sqrt(C[t-1]-R[t]*B[t-1]**2)
  }
  # backward sampling
  theta = matrix(0,nd,n)
  theta[,n] = rnorm(nd,m[n],sqrt(C[n]))
  for (t in (n-1):1)
    theta[,t] = rnorm(nd,m[t]+B[t]*(theta[,t+1]-a[t+1]),H[t])
  if (nd==1){
    theta[1,]
  }
  else{
    theta
  }
}

3.2 Posterior inference

# Prior for beta(0) ~ N(m0,C0)
# Prior for beta(1) ~ N(a1,R1), a1=m0,R1=C0+W
a1     = 0
R1     = 10
av     = 2.01
bv     = 1.01
aw     = 2.01
bw     = 1.01*W    
m0     = 0
C0     = 10-W

# MCMC initual values
V0     = 1
W0     = 0.5
M0     = 10000
M      = 10000

set.seed(65432)
draws1 = mcmc.joint(y,a1,R1,av,bv,aw,bw,M0,M,V0,W0)

# states MCMC autocorrelations and posterior quantiles
q1 = t(apply(draws1[,1:n],2,quantile,c(0.05,0.5,0.95)))
lagmax = 30
acfs = array(0,c(2,n,lagmax))
for (i in 1:n)
  acfs[1,i,] = acf(draws1[,i],plot=FALSE,lag.max=lagmax)$acf[2:(lagmax+1)]

3.2.1 Posterior inference for \(V\) and \(W\)

par(mfrow=c(2,3),mar = c(2, 2, 1, 1))
ts.plot(draws1[,n+1],ylab="",xlab="Iteration",main="V")
acf(draws1[,n+1],main="")
hist(draws1[,n+1],prob=TRUE,main="",xlab="")
ts.plot(draws1[,n+2],ylab="",xlab="Iteration",main="W")
acf(draws1[,n+2],main="")
hist(draws1[,n+2],prob=TRUE,main="",xlab="")

3.2.2 MCMC autocorrelation for states

par(mfrow=c(1,1),mar=c(5, 4, 4, 2))
ts.plot(t(acfs[1,,]),xlab="Lag",ylab="ACF",main="",ylim=range(acfs))

3.2.3 Posterior quantiles for states

par(mfrow=c(1,1))
ts.plot(q1,col=c(3,2,3),ylim=range(y,q1),ylab="States")
points(y)
title("5th, 50th and 95th percentiles")

4 Single-move

The routine mcmc.single samples \(\beta_1,\ldots,\beta_n\), \(V\) and \(W\) iteratively from their univariate full conditional posterior distributions, ie. \[\begin{eqnarray*} p(V|\beta_1,\ldots,\beta_n,W,\mbox{data})\\ p(W|\beta_1,\ldots,\beta_n,V,\mbox{data}), \end{eqnarray*}\] and \[ p(\beta_t|\beta_{-t},V,W,\mbox{data}) \qquad t=1,\ldots,n, \] where \(\beta_{-t}\) is the vector \((\beta_1,\ldots,\beta_n)\) excluding the \(\beta_t\) state.

mcmc.single returns an \(M \times (n+2)\) matrix with \(M\) draws of the \((n+2)\)-dimensional vector \((\beta_1,\ldots,\beta_n,V,W)\).

mcmc.single = function(y,a1,R1,av,bv,aw,bw,M0,M,beta,V,W){
  n         = length(y)
  beta.draw = beta
  V.draw    = V
  W.draw    = W
  niter     = M0+M
  draws     = matrix(0,niter,n+2)
  for (iter in 1:niter){
    beta.draw = single.sampling(y,a1,R1,av,bv,aw,bw,M0,M,beta.draw,V.draw,W.draw)
    par2v     = bv+sum((y-beta.draw)^2)/2
    par2w     = bw+sum((beta.draw[2:n]-beta.draw[1:(n-1)])^2)/2
    V.draw = 1/rgamma(1,av+n/2,par2v)
    W.draw = 1/rgamma(1,aw+(n-1)/2,par2w)
    draws[iter,] = c(beta.draw,V.draw,W.draw)
  }
  return(draws[(M0+1):niter,])
}

4.1 Sampling from \(p(\beta_t|\beta_{-t},V,W,\mbox{data})\)

The routine single.sampling samples, for \(t=1,\ldots,n\), from the univariate full conditional distributons \[ p(\beta_t|\beta_{-t},V,W,\mbox{data}). \]

single.sampling = function(y,a1,R1,av,bv,aw,bw,M0,M,beta,V,W){
  V1 = 1/(1/R1+1/W+1/V)
  b1 = V1*(a1/R1+beta[2]/W+y[1]/V)
  beta[1] = rnorm(1,b1,sqrt(V1))
  for (t in 2:(n-1)){
    Vt = 1/(2/W+1/V)         
    bt = Vt*((beta[t-1]+beta[t+1])/W+y[t]/V)
    beta[t] = rnorm(1,bt,sqrt(Vt))
  }  
  Vn = 1/(1/W+1/V)
  bn = Vn*(beta[n-1]/W+y[n]/V)
  beta[n] = rnorm(1,bn,sqrt(Vn))
  return(beta)
}

4.2 Posterior inference

set.seed(65432)
beta0v = y
draws2 = mcmc.single(y,a1,R1,av,bv,aw,bw,M0,M,beta0v,V0,W0)
q2     = t(apply(draws2[,1:n],2,quantile,c(0.05,0.5,0.95)))
for (i in 1:n)
  acfs[2,i,] = acf(draws2[,i],plot=FALSE,lag.max=lagmax)$acf[2:(lagmax+1)]

4.2.1 Posterior inference for \(V\) and \(W\)

par(mfrow=c(2,3),mar = c(2, 2, 1, 1))
ts.plot(draws2[,n+1],ylab="",xlab="Iteration",main="V")
acf(draws2[,n+1],main="")
hist(draws2[,n+1],prob=TRUE,main="",xlab="")
ts.plot(draws2[,n+2],ylab="",xlab="Iteration",main="W")
acf(draws2[,n+2],main="")
hist(draws2[,n+2],prob=TRUE,main="",xlab="")

4.2.2 MCMC autocorrelation for states

par(mfrow=c(1,1),mar=c(5, 4, 4, 2))
ts.plot(t(acfs[2,,]),xlab="Lag",ylab="ACF",main="",ylim=range(acfs))

4.2.3 Posterior quantiles for states

par(mfrow=c(1,1))
ts.plot(q2,col=c(3,2,3),ylim=range(y,q2),ylab="States")
points(y)
title("5th, 50th and 95th percentiles")

5 Comparing MCMC schemes

5.1 Comparing ACFs

par(mfrow=c(1,1),mar=c(5, 4, 4, 2))
ts.plot(t(acfs[2,,]),xlab="Lag",ylab="ACF",main="",ylim=range(acfs))
for (i in 1:n)
  lines(acfs[1,i,],col=2)
legend("topright",legend=c("Single-move","Block-move"),col=1:2,lwd=2,bty="n")

5.2 Comparing state posterior quantiles

par(mfrow=c(1,1),mar=c(5, 4, 4, 2))
ts.plot(q2)
lines(q1[,1],col=2)
lines(q1[,2],col=2)
lines(q1[,3],col=2)
legend("top",legend=c("Single-move","Block-move"),col=1:2,lwd=2,bty="n")

5.3 Comparing marginal posterior for \(V\) and \(W\)

par(mfrow=c(1,2),mar = c(2, 2, 1, 1))
plot(density(draws2[,n+1]),main="V")
lines(density(draws1[,n+1]),col=2)
plot(density(draws2[,n+2]),main="W")
lines(density(draws1[,n+2]),col=2)
legend("topright",legend=c("Single-move","Block-move"),col=1:2,lwd=2,bty="n")

6 Replications

We replicate the above exercise for 20 sets of time series data and run both MCMC schemes. We assume that for \(n=100\) and $W {0.01,0.25,1.0}.

set.seed(65432123)
nrep  = 20
n     = 100
V     = 1.0
M0    = 1000
M     = 1000
beta0 = 0.0
V0    = 1.0
W0    = 0.5
a1    = 0.0
R1    = 10.0
av    = 2.01
bv    = 1.01
aw    = 2.01
m0    = 0.0
ws    = c(0.01,0.25,1.0)
nws   = length(ws)
neffs = array(0,c(nrep,nws,2))
for (rep in 1:nrep){
  for (i in 1:3){
    y      = sample.data(n,ws[i],V,beta0)
    beta0v = y
    C0     = 10-W
    bw     = 1.01*W    
    draws1 = mcmc.joint(y,a1,R1,av,bv,aw,bw,M0,M,V0,W0)
    acf1  = acf(draws1[,n+1]/draws1[,n+2],lag.max=31,plot=FALSE)
    neff1 = M/(1+2*sum(acf1$acf[2:31]))
    draws2 = mcmc.single(y,a1,R1,av,bv,aw,bw,M0,M,beta0v,V0,W0)
    acf2  = acf(draws2[,n+1]/draws2[,n+2],lag.max=31,plot=FALSE)
    neff2 = M/(1+2*sum(acf2$acf[2:31]))
    neffs[rep,i,] = trunc(c(neff1,neff2))
  }
}

6.1 Effective sample size (ESS)

comb = c("(100,0.01)","(100,0.25)","(100,1)")
par(mfrow=c(1,1),mar=c(5, 4, 4, 2))
boxplot(cbind(neffs[,,1],neffs[,,2]),outline=FALSE,names=c(comb,comb),
        xlab=paste("(n,W) for V=",V,sep=""))
abline(v=3.5)
legend("topleft",legend=c("Block-move"),bty="n")
legend("topright",legend=c("Single-move"),bty="n")
title(paste("MCMC posterior draws = ",M,sep=""))

6.2 Relative ESS

par(mfrow=c(1,1),mar=c(5, 4, 4, 2))
boxplot(neffs[,,2]/neffs[,,1],outline=FALSE,names=comb,ylab="Relative ESS",
        xlab=paste("(n,W) for V=",V,sep=""))
abline(h=1,lty=2)
title("Ratio single-move to block-move")