###################################################################### # # First order DLM - in class R coding # ###################################################################### rm(list=ls()) set.seed(2468) pdf(file="localelvelmodel-inclass.pdf",width=15,height=10) # Simulating a 1st order DLM n = 50 sig2 = 1.0 tau2 = 0.25 x = rep(0,n) y = rep(0,n) y[1] = rnorm(1,x[1],sqrt(sig2)) for (t in 2:n){ x[t] = rnorm(1,x[t-1],sqrt(tau2)) y[t] = rnorm(1,x[t],sqrt(sig2)) } xtrue = x plot(y,xlab="Time",ylab="",pch=16) lines(x,col=2,type="b",pch=16) # Implementing the Kalman filter m0 = 0 C0 = 10 ms = rep(0,n) Cs = rep(0,n) Rs = rep(0,n) As = rep(0,n) Qs = rep(0,n) lvero = rep(0,n) m = m0 C = C0 for (t in 1:n){ Rs[t] = C + tau2 Qs[t] = Rs[t] + sig2 lvero[t] = dnorm(y[t],m,sqrt(Qs[t]),log=TRUE) As[t] = Rs[t]/Qs[t] m = (1-As[t])*m + As[t]*y[t] C = As[t]*sig2 ms[t] = m Cs[t] = C } par(mfrow=c(2,3)) ts.plot(Rs) ts.plot(Qs) ts.plot(lvero) ts.plot(As) ts.plot(ms) ts.plot(Cs) par(mfrow=c(1,1)) L = ms-2*sqrt(Cs) U = ms+2*sqrt(Cs) plot(xtrue,pch=16,ylim=range(xtrue,ms,L,U),xlab="Time",ylab="State") lines(ms,type="b",col=2,pch=16) lines(L,col=2) lines(U,col=2) # Writing the Kalman filter as a function kf = function(y,m0,C0,sig2,tau2){ n = length(y) ms = rep(0,n) Cs = rep(0,n) Rs = rep(0,n) As = rep(0,n) Qs = rep(0,n) lvero = rep(0,n) m = m0 C = C0 for (t in 1:n){ Rs[t] = C + tau2 Qs[t] = Rs[t] + sig2 lvero[t] = dnorm(y[t],m,sqrt(Qs[t]),log=TRUE) As[t] = Rs[t]/Qs[t] m = (1-As[t])*m + As[t]*y[t] C = As[t]*sig2 ms[t] = m Cs[t] = C } return(list(m=ms,C=Cs,L=lvero)) } # Plot of the likelihood function L(sig2,tau2;y) N = 50 sig2s = seq(0.05,3,length=N) tau2s = seq(0.01,2,length=N) lpred = matrix(0,N,N) for (i in 1:N) for (j in 1:N){ runkf = kf(y,m0,C0,sig2s[i],tau2s[j]) lpred[i,j] = sum(runkf$L) } pred = exp(lpred-max(lpred)) contour(sig2s,tau2s,pred,drawlabels=FALSE, xlab=expression(sigma^2), ylab=expression(tau^2)) points(sig2,tau2,pch=16,cex=2,col=2) # Plot of p(x(t)|y(1),...,y(t),sig2,tau2) # --------------------------------------- runkf = kf(y,m0,C0,sig2,tau2) L = runkf$m-2*sqrt(runkf$C) U = runkf$m+2*sqrt(runkf$C) plot(runkf$m,ylim=range(L,U),pch=16,type="b",col=2,xlab="Time",ylab="State") lines(L,col=2) lines(U,col=2) points(xtrue,pch=16) # Running the Kalman smoother # --------------------------- ks = function(y,m0,C0,sig2,tau2){ n = length(y) ms = rep(0,n) Cs = rep(0,n) Rs = rep(0,n) As = rep(0,n) Qs = rep(0,n) lvero = rep(0,n) Bs = rep(0,n) m = m0 C = C0 for (t in 1:n){ Rs[t] = C + tau2 Qs[t] = Rs[t] + sig2 lvero[t] = dnorm(y[t],m,sqrt(Qs[t]),log=TRUE) As[t] = Rs[t]/Qs[t] m = (1-As[t])*m + As[t]*y[t] C = As[t]*sig2 ms[t] = m Cs[t] = C Bs[t] = C/(C+tau2) } mn = rep(0,n) Cn = rep(0,n) mn[n] = ms[n] Cn[n] = Cs[n] for (t in (n-1):1){ mn[t] = (1-Bs[t])*ms[t]+Bs[t]*mn[t+1] Cn[t] = (1-Bs[t])*Cs[t]+Bs[t]^2*Cn[t+1] } return(list(m=ms,C=Cs,L=lvero,mb=mn,Cb=Cn)) } # Plot of p(x(t)|y(1),...,y(n),sig2,tau2) runks = ks(y,m0,C0,sig2,tau2) L = runks$mb-2*sqrt(runks$Cb) U = runks$mb+2*sqrt(runks$Cb) Lf = runks$m-2*sqrt(runks$C) Uf = runks$m+2*sqrt(runks$C) plot(runks$mb,ylim=range(L,U,Lf,Uf), pch=16,type="b",col=2,ylab="",xlab="Time") lines(L,col=2) lines(U,col=2) lines(Lf,col=4) lines(Uf,col=4) lines(runks$m,col=4,type="b",pch=16) # Implementing the FFBS ffbs = function(y,m0,C0,sig2,tau2){ n=length(y);ms=rep(0,n);Bs=rep(0,n);m=m0;C=C0 for (t in 1:n){ R = C + tau2;Q = R + sig2 A = R/Q m = (1-A)*m + A*y[t] C = A*sig2 ms[t] = m; Bs[t] = C/(C+tau2) } a = ms[n];R = C x = rep(0,n) x[n] = rnorm(1,a,sqrt(R)) for (t in (n-1):1){ a = (1-Bs[t])*ms[t]+Bs[t]*x[t+1] R = Bs[t]*tau2 x[t] = rnorm(1,a,sqrt(R)) } return(x) } # Testing the FFBS function xdraw = matrix(0,1000,n) for (i in 1:1000) xdraw[i,] = ffbs(y,m0,C0,sig2,tau2) qx = apply(xdraw,2,quantile,c(0.025,0.5,0.975)) ts.plot(t(qx)) lines(L,col=2) lines(U,col=2) lines(runks$mb,col=2) # THE FULL MCMC (GIBBS SAMPLER) IN ACTION set.seed(5555) a0 = 0.01;b0=0.01;c0=0.01;d0=0.01 niter = 1000 xs = matrix(0,niter,n) sig2s = rep(0,niter) tau2s = rep(0,niter) sig2 = 5 tau2 = 3 system.time({ for (iter in 1:niter){ # sample the states x1,...,xn x = ffbs(y,m0,C0,sig2,tau2) # sample sig2 par1 = a0 + n/2 par2 = b0 + sum((y-x)^2)/2 sig2 = 1/rgamma(1,par1,par2) # sample tau2 par1 = c0 + n/2 par2 = d0 + sum(diff(x)^2)/2 tau2 = 1/rgamma(1,par1,par2) # store the samples xs[iter,] = x sig2s[iter] = sig2 tau2s[iter] = tau2 } }) xs.ffbs = xs sig2s.ffbs = sig2s tau2s.ffbs = tau2s qx = apply(xs,2,quantile,c(0.025,0.5,0.975)) par(mfrow=c(2,3)) ts.plot(t(qx),ylab="State") acf(sig2s,main=expression(sigma^2)) hist(sig2s,main=expression(sigma^2),prob=TRUE,xlab="") acf(tau2s,main=expression(tau^2)) hist(tau2s,main=expression(tau^2),prob=TRUE,xlab="") plot(sig2s,tau2s,xlab=expression(sigma^2),ylab=expression(tau^2),pch=16) par(mfrow=c(1,1)) acf1 = acf(xs[,1],plot=FALSE) plot(acf1$acf,type="l",ylab="ACF",xlab="LAG") for (i in 2:n){ acf1 = acf(xs[,i],plot=FALSE) lines(acf1$acf) } # THE FULL MCMC - SINGLE MOVE set.seed(5555) a0 = 0.01;b0=0.01;c0=0.01;d0=0.01 niter = 10000 xs = matrix(0,niter,n) sig2s = rep(0,niter) tau2s = rep(0,niter) sig2 = 5 tau2 = 3 x = xtrue system.time({ for (iter in 1:niter){ # sample the state x1 var = 1/(1/tau2+1/sig2) mean = var*(x[2]/tau2+y[1]/sig2) x[1] = rnorm(1,mean,sqrt(var)) # sample the state xn var = 1/(1/tau2+1/sig2) mean = var*(x[n-1]/tau2+y[n]/sig2) x[n] = rnorm(1,mean,sqrt(var)) # sample the states x2,....,xn-1 for (t in 2:(n-1)){ var = 1/(2/tau2+1/sig2) mean = var*((x[t-1]+x[t+1])/tau2+y[t]/sig2) x[t] = rnorm(1,mean,sqrt(var)) } # sample sig2 par1 = a0 + n/2 par2 = b0 + sum((y-x)^2)/2 sig2 = 1/rgamma(1,par1,par2) # sample tau2 par1 = c0 + n/2 par2 = d0 + sum(diff(x)^2)/2 tau2 = 1/rgamma(1,par1,par2) # store the samples xs[iter,] = x sig2s[iter] = sig2 tau2s[iter] = tau2 } }) xs.single = xs sig2s.single = sig2s tau2s.single = tau2s par(mfrow=c(1,2)) acf1 = acf(xs.ffbs[,1],plot=FALSE) plot(acf1$acf,type="l",ylim=c(-0.1,1),xlab="LAG",ylab="ACF of states") for (i in 2:n){ acf1 = acf(xs.ffbs[,i],plot=FALSE) lines(acf1$acf) } title("Sampling x1,...,xn jointly") acf1 = acf(xs.single[,1],plot=FALSE) plot(acf1$acf,type="l",ylim=c(-0.1,1),xlab="LAG",ylab="ACF of states") for (i in 2:n){ acf1 = acf(xs.single[,i],plot=FALSE) lines(acf1$acf) } title("Sampling x1,...,xn individually") dev.off()