############################################################################################ # # FIRST ORDER DYNAMIC LINEAR MODEL # ############################################################################################ # # 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/ # ############################################################################################ rm(list=ls()) DLM = function(y,mu,beta,sig2,tau2,m0,C0){ n = length(y) m = m0 C = C0 l = 0.0 for (t in 1:n){ a = mu+beta*m R = beta^2*C+tau2 Q = R+sig2 l = l+dnorm(y[t],a,sqrt(Q),log=TRUE) A = R/Q m = (1-A)*a+A*y[t] C = A*sig2 } return(l) } inv = function(s){ A = matrix(c(s[1],s[2],s[2],s[3]),2,2) iA = solve(A) draw = t(chol(iA))%*%rnorm(2) return(c(iA[1,1],iA[1,2],iA[2,2],draw[1],draw[2])) } prod = function(s){ A = matrix(c(s[1],s[2],s[2],s[3]),2,2) x = matrix(c(s[4],s[5]),2,1) b = A%*%x return(c(b[1],b[2])) } LW = function(y,sig,tau,delta,b0,sB0,m0,sC0,N){ n = length(y) h2 = 1-((3*delta-1)/(2*delta))^2 a = sqrt(1-h2) pars = matrix(b0+sB0*rnorm(2*N),N,2) xs = rnorm(N,m0,sC0) xss = NULL parss = array(0,c(N,2,n)) for (t in 1:n){ mpar = apply(pars,2,mean) vpar = var(pars) ms = a*pars+(1-a)*matrix(mpar,N,2,byrow=T) mus = pars[,1]+pars[,2]*xs w = dnorm(y[t],mus,sig) k = sample(1:N,size=N,replace=T,prob=w) ms1 = ms[k,] + matrix(rnorm(2*N),N,2)%*%chol(h2*vpar) xt = ms[,1] + ms[,2]*xs[k] + tau*rnorm(N) w1 = dnorm(y[t],xt,sig)/w[k] k = sample(1:N,size=N,replace=T,prob=w1) xs = xt[k] pars = ms1[k,] xss = rbind(xss,xs) parss[,,t] = pars } return(list(pars=parss,xs=xss)) } LWFA = function(y,sig,tau,delta,b0,sB0,m0,sC0,N){ n = length(y) h2 = 1-((3*delta-1)/(2*delta))^2 a = sqrt(1-h2) sd = sqrt(sig2+tau2) var = 1/(1/sig2+1/tau2) sd1 = sqrt(var) pars = matrix(b0+sB0*rnorm(2*N),N,2) xs = rnorm(N,m0,sC0) xss = NULL parss = array(0,c(N,2,n)) for (t in 1:n){ mpar = apply(pars,2,mean) vpar = var(pars) ms = a*pars+(1-a)*matrix(mpar,N,2,byrow=T) mus = pars[,1]+pars[,2]*xs w = dnorm(y[t],mus,sd) k = sample(1:N,size=N,replace=T,prob=w) pars = ms[k,] + matrix(rnorm(2*N),N,2)%*%chol(h2*vpar) mean = var*(y[t]/sig2+(pars[,1]+pars[,2]*xs)/tau2) xs = mean + sd1*rnorm(N) xss = rbind(xss,xs) parss[,,t] = pars } return(list(pars=parss,xs=xss)) } APFSS= function(y,sig,tau,b0,B0,m0,C0,N){ n = length(y) sB0 = sqrt(B0) sC0 = sqrt(C0) pars = matrix(b0+sB0*rnorm(2*N),N,2) xs = rnorm(N,m0,sC0) xss = NULL parss = array(0,c(N,2,n)) s = matrix(0,N,5) s[,1] = 1/C0 s[,2] = 0.0 s[,3] = 1/C0 s[,4] = b0/C0 s[,5] = b0/C0 for (t in 1:n){ w = dnorm(y[t],pars[,1]+pars[,2]*xs,sig) k = sample(1:N,size=N,replace=T,prob=w) pars = pars[k,] xs = xs[k] xt = rnorm(N,pars[,1]+pars[,2]*xs,tau) w1 = dnorm(y[t],xt,sig)/w[k] k = sample(1:N,size=N,replace=T,prob=w1) xs1 = xt[k] s[,1] = s[k,1] + 1 s[,2] = s[k,2] + xs s[,3] = s[k,3] + xs^2 s[,4] = s[k,4] + xs1 s[,5] = s[k,5] + xs1*xs xs = xs1 mat = cbind(s[,1],s[,2],s[,3]) mat1 = t(apply(mat,1,inv)) draws = mat1[,4:5] mat = cbind(mat1[,1:3],s[,4],s[,5]) coef = t(apply(mat,1,prod)) pars = coef + draws xss = rbind(xss,xs) parss[,,t] = pars } return(list(pars=parss,xs=xss)) } PL = function(y,sig2,tau2,b0,B0,m0,C0,N){ n = length(y) sB0 = sqrt(B0) sC0 = sqrt(C0) sd = sqrt(sig2+tau2) var = 1/(1/sig2+1/tau2) sd1 = sqrt(var) pars = matrix(b0+sB0*rnorm(2*N),N,2) xs = rnorm(N,m0,sC0) xss = NULL parss = array(0,c(N,2,n)) s = matrix(0,N,5) s[,1] = 1/C0 s[,2] = 0.0 s[,3] = 1/C0 s[,4] = b0/C0 s[,5] = b0/C0 for (t in 1:n){ w = dnorm(y[t],pars[,1]+pars[,2]*xs,sd) k = sample(1:N,size=N,replace=T,prob=w) xs = xs[k] pars = pars[k,] mean = var*(y[t]/sig2+(pars[,1]+pars[,2]*xs)/tau2) xs1 = mean + sd1*rnorm(N) s[,1] = s[k,1] + 1 s[,2] = s[k,2] + xs s[,3] = s[k,3] + xs^2 s[,4] = s[k,4] + xs1 s[,5] = s[k,5] + xs1*xs xs = xs1 mat = cbind(s[,1],s[,2],s[,3]) mat1 = t(apply(mat,1,inv)) draws = mat1[,4:5] mat = cbind(mat1[,1:3],s[,4],s[,5]) coef = t(apply(mat,1,prod)) pars = coef + draws xss = rbind(xss,xs) parss[,,t] = pars } return(list(pars=parss,xs=xss)) } # Simulating 1st order DLM # ------------------------ set.seed(2468) n = 10 x0 = 0.0 mu = 1.0 beta = 0.95 sig2 = 1.00 tau2 = 0.50 sig = sqrt(sig2) tau = sqrt(tau2) x = rep(0,n) x[1] = rnorm(1,mu+beta*x0,tau) for (t in 2:n) x[t] = rnorm(1,mu+beta*x[t-1],tau) y = rnorm(n,x,sig) pdf(file="lw-lwfa-apfss-pl.pdf",width=14,height=8) par(mfrow=c(1,1)) plot(y,pch=16,xlab="Time",ylab="",main="y(t) and x(t)",cex=1) lines(x,col=2,lwd=3,type="b",pch=16) # Prior information # ----------------- b0 = 0.0 B0 = 1 m0 = 0.0 C0 = 1 sB0 = sqrt(B0) sC0 = sqrt(C0) # True posterior distribution p(beta|data) # ---------------------------------------- ts = c(1,2,3,8,9,10) nt = length(ts) N = 100 mus = seq(-2,2.5,length=N) betas = seq(-2,2.5,length=N) hbeta = betas[2]-betas[1] post = array(0,c(nt,N,N)) for (i in 1:N){ for (j in 1:N){ post[,i,j] = dnorm(mus[i],b0,sB0,log=TRUE)+dnorm(betas[j],b0,sB0,log=TRUE) for (t in 1:nt) post[t,i,j] = post[t,i,j] + DLM(y[1:ts[t]],mus[i],betas[j],sig2,tau2,m0,C0) } } post = exp(post-max(post)) for (t in 1:nt) post[t,,] = post[t,,]/max(post[t,,]) tt = c(1,nt) mpost = array(0,c(2,2,N)) for (i in 1:2) for (j in 1:2){ mpost[i,j,] = apply(post[tt[i],,],j,sum) mpost[i,j,] = mpost[i,j,]/max(mpost[i,j,]) } par(mfrow=c(2,3)) for (t in 1:nt){ contour(mus,betas,post[t,,],xlab=expression(mu),ylab=expression(beta),drawlabels=FALSE) title(paste("time t=",ts[t],sep="")) } names = c("mu","beta") parss = cbind(mus,betas) par(mfrow=c(2,2)) for (i in 1:2) for (j in 1:2){ plot(parss[,i],mpost[j,i,],xlab="",ylab="",main="",type="l") title(paste("p(",names[i],"|y^",ts[tt[j]],")",sep="")) } # Sequential Monte Carlo Methods # ------------------------------ N = 1000 deltas = c(0.75,0.85,0.95,0.99) L1 = min(mus) U1 = max(mus) L2 = min(betas) U2 = max(betas) # Liu-West filter # --------------- set.seed(13579) for (delta in deltas){ lw = LW(y,sig,tau,delta,b0,sB0,m0,sC0,N) par(mfrow=c(2,3)) for (t in 1:nt){ plot(lw$pars[,,ts[t]],col=2,pch=16,cex=0.5,xlab=expression(mu),ylab=expression(beta), xlim=c(L1,U1),ylim=c(L2,U2)) contour(mus,betas,post[t,,], drawlabels=FALSE,add=TRUE) title(paste("LW: time t=",ts[t],"\n delta=",delta,sep="")) } } # Fully adapted Liu-West # ---------------------- set.seed(13579) for (delta in deltas){ lwfa = LWFA(y,sig,tau,delta,b0,sB0,m0,sC0,N) par(mfrow=c(2,3)) for (t in 1:nt){ plot(lwfa$pars[,,ts[t]],col=2,pch=16,cex=0.5,xlab=expression(mu),ylab=expression(beta), xlim=c(L1,U1),ylim=c(L2,U2)) contour(mus,betas,post[t,,], drawlabels=FALSE,add=TRUE) title(paste("LWFA: time t=",ts[t],"\n delta=",delta,sep="")) } } # Auxiliary particle filter + sufficient statistics # ------------------------------------------------- set.seed(13579) for (delta in deltas){ apfss = APFSS(y,sig,tau,b0,B0,m0,C0,N) par(mfrow=c(2,3)) for (t in 1:nt){ plot(apfss$pars[,,ts[t]],col=2,pch=16,cex=0.5,xlab=expression(mu),ylab=expression(beta), xlim=c(L1,U1),ylim=c(L2,U2)) contour(mus,betas,post[t,,], drawlabels=FALSE,add=TRUE) title(paste("APFSS: time t=",ts[t],"\n delta=",delta,sep="")) } } # Particle Learning # ----------------- set.seed(13579) pl = PL(y,sig2,tau2,b0,B0,m0,C0,N) par(mfrow=c(2,3)) for (t in 1:nt){ plot(pl$pars[,,ts[t]],col=2,pch=16,cex=0.5,xlab=expression(mu),ylab=expression(beta), xlim=c(L1,U1),ylim=c(L2,U2)) contour(mus,betas,post[t,,], drawlabels=FALSE,add=TRUE) title(paste("PL: time t=",ts[t],"\n PL",sep="")) } names = c("mu","beta") parss = cbind(mus,betas) j = 2 par(mfrow=c(1,1)) for (i in 1:2){ plot(parss[,i],mpost[j,i,],xlab="",ylab="",main="",type="l",lwd=2) title(paste("p(",names[i],"|y^",ts[tt[j]],")",sep="")) dlw = density(lw$pars[,i,ts[nt]]) lines(dlw$x,dlw$y/max(dlw$y),col=2,lwd=2) dlwfa = density(lwfa$pars[,i,ts[nt]]) lines(dlwfa$x,dlwfa$y/max(dlwfa$y),col=3,lwd=2) dapfss = density(apfss$pars[,i,ts[nt]]) lines(dapfss$x,dapfss$y/max(dapfss$y),col=4,lwd=2) dpl = density(pl$pars[,i,ts[nt]]) lines(dpl$x,dpl$y/max(dpl$y),col=6,lwd=2) lines(parss[,i],mpost[j,i,],lwd=2,col=1) legend(-2,1,legend=c("TRUE","LW","LWFA","APFSS","PL"),col=c(1:4,6), lty=rep(1,5),lwd=rep(2,5),bty="n") } dev.off()