############################################################################################ # # 1st ORDER DYNAMIC LINEAR MODEL # # For t=1,...,n, # # y(t) ~ N(theta(t);V) # theta(t) ~ N(theta(t-1);W) # # Prior distributions: # # theta(0) ~ N(m0,C0) # V ~ IG(a,b) # W ~ IG(c,d) # ############################################################################################ # # HEDIBERT FREITAS LOPES # Associate Professor of Econometrics and Statistics # The University of Chicago Booth School of Business # 5807 South Woodlawn Avenue # Chicago, Illinois, 60637 # Email : hlopes@ChicagoGSB.edu # URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/ # ############################################################################################ ffbs = function(y,V,W,m0,C0){ n = length(y) 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) mm = rep(0,n) CC = rep(0,n) x = rep(0,n) llike = 0.0 for (t in 1:n){ if(t==1){ a[1] = m0 R[1] = C0 + W }else{ a[t] = m[t-1] R[t] = C[t-1] + W } f = a[t] Q = R[t] + V A = R[t]/Q m[t] = a[t]+A*(y[t]-f) C[t] = R[t]-Q*A**2 B[t-1] = C[t-1]/R[t] H[t-1] = C[t-1]-R[t]*B[t-1]**2 llike = llike + dnorm(y[t],f,sqrt(Q),log=TRUE) } mm[n] = m[n] CC[n] = C[n] x[n] = rnorm(1,m[n],sqrt(C[n])) for (t in (n-1):1){ mm[t] = m[t] + C[t]/R[t+1]*(mm[t+1]-a[t+1]) CC[t] = C[t] - (C[t]^2)/(R[t+1]^2)*(R[t+1]-CC[t+1]) x[t] = rnorm(1,m[t]+B[t]*(x[t+1]-a[t+1]),sqrt(H[t])) } return(list(x=x,m=m,C=C,mm=mm,CC=CC,llike=llike)) } q025 = function(x){quantile(x,0.025)} q975 = function(x){quantile(x,0.975)} # Simulating 1st order DLM # ------------------------ set.seed(123456) W = 0.5 V = 1.0 n = 100 m0 = 0.0 C0 = 10.0 x0 = 0.00 w = rnorm(n,0,sqrt(W)) v = rnorm(n,0,sqrt(V)) x = y = rep(0,n) x[1] = x0 + w[1] y[1] = x[1] + v[1] for (t in 2:n){ x[t] = x[t-1] + w[t] y[t] = x[t] + v[t] } run = ffbs(y,V,W,m0,C0) m = run$m C = run$C mm = run$mm CC = run$CC L1 = m-2*C U1 = m+2*C L2 = mm-2*CC U2 = mm+2*CC pdf(file="local-level-sim.pdf",width=10,height=7) par(mfrow=c(1,1)) plot(y1,type="l",xlab="time",ylab="",lwd=3,ylim=c(min(L1,L2,y,x),max(U1,U2,y,x))) lines(x1,col=2,lwd=3) legend(70,8,legend=c("y(t)","x(t)"),lty=c(1,1),col=c(1,2),lwd=c(3,3),bty="n") dev.off() pdf(file="local-level-ffbs.pdf",width=10,height=7) par(mfrow=c(1,1)) ts.plot(m,xlab="time",ylab="",ylim=c(min(L1,L2,y,x),max(U1,U2,y,x)),lwd=2) lines(L1,lwd=2) lines(U1,lwd=2) lines(mm,col=2,lwd=2) lines(L2,col=2,lwd=2) lines(U2,col=2,lwd=2) legend(70,8,legend=c("Forward filtering","Backward smoothing"), lty=c(1,1),col=c(1,2),lwd=c(3,3),bty="n") dev.off() N = 50 Vs = seq(0.1,2,length=N) Ws = seq(0.1,2,length=N) likes = matrix(0,N,N) for (i in 1:N){ for (j in 1:N){ V = Vs[i] W = Ws[j] run = ffbs(y,V,W,m0,C0) likes[i,j] = run$llike } } pdf(file="local-level-contour.pdf",width=10,height=7) contour(Vs,Ws,exp(likes),xlab=expression(sigma^2),ylab=expression(tau^2),drawlabels=FALSE) dev.off() # Hyperparameters a = 0.01 b = 0.01 c = 0.01 d = 0.01 # MCMC step set.seed(1255) burn = 1000 M = 1000 niter = burn + M V1 = V W1 = W draws = NULL for (iter in 1:niter){ run = ffbs(y,V1,W1,m0,C0) x = run$x V1 = 1/rgamma(1,a+n/2,b+sum((y-x)^2)/2) W1 = 1/rgamma(1,c+(n-1)/2,d+sum(diff(x)^2)/2) draws = rbind(draws,c(V1,W1,x)) } draws = draws[(burn+1):(niter),] xs = draws[,3:(n+2)] lx = apply(xs,2,q025) mx = apply(xs,2,median) ux = apply(xs,2,q975) pdf(file="sig2tau2.pdf",width=13,height=8) par(mfrow=c(1,3)) contour(Vs,Ws,exp(likes),xlab=expression(sigma^2),ylab=expression(tau^2),drawlabels=FALSE) points(draws[,1:2],pch=16,col=2,cex=1) hist(draws[,1],xlab="",ylab="Density",main=expression(sigma^2)) hist(draws[,2],xlab="",ylab="Density",main=expression(tau^2)) dev.off() pdf(file="local-level-ffbsx.pdf",width=10,height=7) par(mfrow=c(1,1)) ts.plot(cbind(lx,mx,ux),xlab="time",ylab="",main="") dev.off()