################################################################################################ # # Simple normal dynamic linear regression # # Block sampling the states vs single sampling the states # # By Hedibert F. Lopes # March 14th 2021 # ################################################################################################ sim.simplendlreg = function(n,mu,phi,tau2,sig2,beta0){ tau = sqrt(tau2) sig = sqrt(sig2) beta = rep(0,n) beta[1] = rnorm(1,mu+phi*beta0,tau) for (t in 2:n) beta[t] = rnorm(1,mu+phi*beta[t-1],tau) x = rnorm(n) y = rnorm(n,beta*x,sig) return(list(x=x,y=y,beta=beta)) } ffbs = function(y,x,m0,C0,mu,phi,tau2,sig2){ n = length(y) m = m0 C = C0 phi2 = phi^2 x2 = x^2 a = rep(0,n) R = rep(0,n) mf = rep(0,n) Cf = rep(0,n) for (t in 1:n){ a[t] = mu + phi*m R[t] = phi2*C + tau2 f = x[t]*a[t] Q = x2[t]*R[t] + sig2 A = R[t]*x[t]/Q m = a[t]+A*(y[t]-f) C = R[t] - A^2*Q mf[t] = m Cf[t] = C } beta = rep(0,n) beta[n] = rnorm(1,mf[n],Cf[n]) for (t in (n-1):1){ var = 1/(phi2/tau2+1/Cf[t]) mean = var*(beta[t+1]*phi/tau2+mf[t]/Cf[t]) beta[t] = rnorm(1,mean,sqrt(var)) } return(beta) } singlemove.beta = function(y,x,m0,C0,mu,phi,tau2,sig2,beta){ b0 = mu+phi*m0 V0 = tau2+phi^2*C0 var = 1/(1/V0+phi^2/tau2+x[1]^2/sig2) mean = var*(b0/V0+phi*(beta[2]-mu)/tau2+x[1]*y[1]/sig2) beta[1] = rnorm(1,mean,sqrt(var)) for (t in 2:(n-1)){ var = 1/((1+phi^2)/tau2+x[t]^2/sig2) mean = var*((mu*(1-phi)+phi*(beta[t-1]+beta[t+1])/tau2)+y[t]*x[t]/sig2) beta[t] = rnorm(1,mean,sqrt(var)) } var = 1/(1/tau2+x[n]^2/sig2) mean = var*((mu+phi*beta[n-1])/tau2+x[n]*y[b=n]/sig2) beta[n] = rnorm(1,mean,sqrt(var)) return(beta) } blocksampling.MCMC = function(y,x,mu,phi,tau2,sig2,m0,C0,burnin,ndraw,thin){ n = length(y) niter = burnin + thin*ndraw pars = matrix(0,niter,4) betas = matrix(0,niter,n) for (iter in 1:niter){ beta = ffbs(y,x,m0,C0,mu,phi,tau2,sig2) tau2 = 1/rgamma(1,(n-1)/2,sum((beta[2:n]-mu-phi*beta[1:(n-1)])^2)/2) sig2 = 1/rgamma(1,n/2,sum((y-x*beta)^2)/2) y1 = beta[2:n]-mu x1 = beta[1:(n-1)] var = tau2/sum(x1^2) mean = sum(x1*y1)/sum(x1^2) phi = rnorm(1,mean,sqrt(var)) y1 = beta[2:n]-phi*beta[1:(n-1)] mu = rnorm(1,mean(y1),sqrt(tau2/(n-1))) pars[iter,] = c(mu,phi,tau2,sig2) betas[iter,] = beta } ind = seq(burnin+1,niter,by=thin) pars = pars[ind,] betas = betas[ind,] return(list(pars=pars,betas=betas)) } singlesampling.MCMC = function(y,x,mu,phi,tau2,sig2,beta,m0,C0,burnin,ndraw,thin){ n = length(y) niter = burnin + thin*ndraw pars = matrix(0,niter,4) betas = matrix(0,niter,n) for (iter in 1:niter){ beta = singlemove.beta(y,x,m0,C0,mu,phi,tau2,sig2,beta) tau2 = 1/rgamma(1,(n-1)/2,sum((beta[2:n]-mu-phi*beta[1:(n-1)])^2)/2) sig2 = 1/rgamma(1,n/2,sum((y-x*beta)^2)/2) y1 = beta[2:n]-mu x1 = beta[1:(n-1)] var = tau2/sum(x1^2) mean = sum(x1*y1)/sum(x1^2) phi = rnorm(1,mean,sqrt(var)) y1 = beta[2:n]-phi*beta[1:(n-1)] mu = rnorm(1,mean(y1),sqrt(tau2/(n-1))) pars[iter,] = c(mu,phi,tau2,sig2) betas[iter,] = beta } ind = seq(burnin+1,niter,by=thin) pars = pars[ind,] betas = betas[ind,] return(list(pars=pars,betas=betas)) } ################################################################################################# # Simulating the data set.seed(31416) n = 100 sig = 1.0 tau = 0.2 mu = 0.0 phi = 1.0 beta0 = 1.5 tau2 = tau^2 sig2 = sig^2 sim = sim.simplendlreg(n,mu,phi,tau2,sig2,beta0) x = sim$x y = sim$y beta = sim$beta param = c(mu,phi,tau2,sig2) pdf(file="NDLM-blocksampling-vs-singlesampling.pdf",width=12,height=9) par(mfrow=c(1,2)) plot(beta,main="Time-varying coefficient",xlab="Time",ylab=expression(beta),type="l") plot(x,y,xlim=range(x),ylim=range(y),xlab="X",ylab="Y") # MCMC set-up set.seed(3251241) m0 = 0 C0 = 10 burnin = 10000 ndraw = 10000 thin = 1 # Joint sampling beta[1],...,beta[n] given (y[1],x[1]),...,(y[n],x[n]) runtime = system.time({ run.mcmc = blocksampling.MCMC(y,x,mu,phi,tau2,sig2,m0,C0,burnin,ndraw,thin) }) pars = run.mcmc$pars betas = run.mcmc$betas qbeta = t(apply(betas,2,quantile,c(0.025,0.5,0.975))) # Individual sampling beta[1],...,beta[n] given (y[1],x[1]),...,(y[n],x[n]) beta = qbeta[,2] runtime1 = system.time({ run.mcmc1 = singlesampling.MCMC(y,x,mu,phi,tau2,sig2,beta,m0,C0,burnin,ndraw,thin) }) pars1 = run.mcmc1$pars betas1 = run.mcmc1$betas qbeta1 = t(apply(betas1,2,quantile,c(0.025,0.5,0.975))) # Comparison names = c(expression(mu),expression(phi),expression(tau^2),expression(sigma^2)) par(mfrow=c(3,4)) for (i in 1:4){ ts.plot(pars[,i],xlab="Iterations",ylab="",main=names[i]) lines(pars1[,i],col=2) } L = 50 for (i in 1:4){ acf = acf(pars[,i],plot=F,lag.max=L) acf1 = acf(pars1[,i],plot=F,lag.max=L) range = range(acf$acf[2:(L+1)],acf1$acf[2:(L+1)]) plot(acf$lag[2:(L+1)],acf$acf[2:(L+1)],xlab="Lag",ylab="ACF",type="l",ylim=range,axes=FALSE) axis(2);box();axis(1,at=seq(1,L,by=3),lab=seq(1,L,by=3)) lines(acf1$lag[2:(L+1)],acf1$acf[2:(L+1)],col=2) abline(h=0,lty=2) } for (i in 1:4){ den = density(pars[,i]) den1 = density(pars1[,i]) plot(den$x,den$y,ylim=range(den$y,den1$y),type="l",ylab="Density",xlab="") lines(den1$x,den1$y,col=2) } ts = trunc(seq(1,n,length=9)) par(mfrow=c(3,3)) for (t in ts){ den = density(betas[,t]) den1 = density(betas1[,t]) plot(den$x,den$y,ylim=range(den$y,den1$y),type="l",ylab="Density",xlab="") lines(den1$x,den1$y,col=2) title(paste("beta(",t,")",sep="")) } L = 100 par(mfrow=c(1,1)) acf = acf(betas1[,1],plot=FALSE,lag.max=L) plot(acf$lag[2:(L+1)],acf$acf[2:(L+1)],xlab="Lag",ylab="ACF",type="l",axes=FALSE,ylim=c(-0.1,1.0)) axis(2);box();axis(1,at=seq(1,L,by=3),lab=seq(1,L,by=3)) for (i in 2:n){ acf = acf(betas1[,i],plot=FALSE,lag.max=L) lines(acf$lag[2:(L+1)],acf$acf[2:(L+1)]) } for (i in 1:n){ acf = acf(betas[,i],plot=FALSE,lag.max=L) lines(acf$lag[2:(L+1)],acf$acf[2:(L+1)],col=2) } legend("topright",legend=c("Block sampling (FFBS)","Individual sampling"),col=2:1,lty=1,bty="n",lwd=2) dev.off()