############################################################################################ # # SEQUENTIAL MONTE CARLO METHODS # # 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/ # ############################################################################################ DLM = function(y,sig2,tau2,m0,C0){ n = length(y) m = rep(0,n) C = rep(0,n) for (t in 1:n){ if (t==1){ a = m0 R = C0 + tau2 }else{ a = m[t-1] R = C[t-1] + tau2 } A = R/(R+sig2) m[t] = (1-A)*a + A*y[t] C[t] = A*sig2 } return(list(m=m,C=C)) } dlm.sim = function(n,sig,tau,x0){ x = rep(0,n) y = rep(0,n) x[1] = rnorm(1,x0,tau) y[1] = rnorm(1,x[1],sig) for (t in 2:n){ x[t] = rnorm(1,x[t-1],tau) y[t] = rnorm(1,x[t],sig) } return(list(x=x,y=y)) } q025=function(x){quantile(x,0.025)} q975=function(x){quantile(x,0.975)} # Simulating 1st order DLM # ------------------------ set.seed(123456) n = 50 tau2 = 0.5 sig2 = c(0.5,1,2)*tau2 nsig = length(sig2) x0 = 0.0 sig = sqrt(sig2) tau = sqrt(tau2) x = matrix(0,n,nsig) y = matrix(0,n,nsig) for (i in 1:3){ sim = dlm.sim(n,sig[i],tau,x0) x[,i] = sim$x y[,i] = sim$y } # Prior hyperparameters # --------------------- m0 = 0 C0 = 100 # Closed form online estimation # ----------------------------- m = matrix(0,n,nsig) C = matrix(0,n,nsig) L = matrix(0,n,nsig) U = matrix(0,n,nsig) for (i in 1:nsig){ dlm = DLM(y[,i],sig2[i],tau2,m0,C0) m[,i] = dlm$m C[,i] = dlm$C L[,i] = m[,i] + qnorm(0.025)*sqrt(C[,i]) U[,i] = m[,i] + qnorm(0.975)*sqrt(C[,i]) } pdf(file="online.pdf",width=15,height=10) ind = c(1,seq(10,n,by=10)) par(mfrow=c(2,nsig)) for (i in 1:nsig){ L1 = min(L[,i],y[,i],x[,i]) U1 = max(U[,i],y[,i],x[,i]) plot(m[,i],xlab="time",ylab="",ylim=c(L1,U1),lwd=2,axes=F,pch=16) axis(2);axis(1,at=ind,label=ind);box() lines(x[,i],lwd=2,col=2) title(paste("tau/sig=",round(sqrt(tau2/sig2[i]),3),sep="")) } for (i in 1:nsig){ L1 = min(L[,i],y[,i],x[,i]) U1 = max(U[,i],y[,i],x[,i]) plot(m[,i],xlab="time",ylab="",ylim=c(L1,U1),lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() lines(L[,i],lwd=2) lines(U[,i],lwd=2) } dev.off() # Sequential importance sampling (SIS) # ------------------------------------ N = 1000 wss = array(0,c(nsig,n,N)) xss = array(0,c(nsig,n,N)) Ess = matrix(0,nsig,n) for (i in 1:nsig){ set.seed(12545) x = rnorm(N,m0,sqrt(C0)) w = rep(1/N,N) ws = NULL xs = NULL ess = NULL for (t in 1:n){ x = rnorm(N,x,tau) w = w*dnorm(y[t,i],x,sig[i]) w1 = w/sum(w) ESS = 1/sum(w1^2) draw = sample(x,size=N,replace=T,prob=w1) xs = rbind(xs,draw) ws = rbind(ws,w1) ess = c(ess,ESS) } wss[i,,] = ws xss[i,,] = xs Ess[i,] = ess } mx = matrix(0,n,nsig) Lx = matrix(0,n,nsig) Ux = matrix(0,n,nsig) for (i in 1:nsig){ mx[,i] = apply(xss[i,,],1,median) Lx[,i] = apply(xss[i,,],1,q025) Ux[,i] = apply(xss[i,,],1,q975) } pdf(file=paste("sis",N,".pdf",sep=""),width=15,height=10) ind = c(1,seq(10,n,by=10)) par(mfrow=c(2,nsig)) for (i in 1:nsig){ L1 = min(L[,i],Lx[,i]) U1 = max(U[,i],Ux[,i]) plot(m[,i],xlab="time",ylab="",ylim=c(L1,U1),lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() lines(L[,i],lwd=2) lines(U[,i],lwd=2) lines(Lx[,i],lwd=2,col=2) lines(mx[,i],lwd=2,col=2) lines(Ux[,i],lwd=2,col=2) } for (i in 1:nsig){ plot(Ess[i,],xlab="time",ylab="",lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() } dev.off() # sequential importance sampling with resampling (SISR) or Bayesian bootstrap filter # ---------------------------------------------------------------------------------- N = 1000 wss = array(0,c(nsig,n,N)) xss = array(0,c(nsig,n,N)) Ess = matrix(0,nsig,n) for (i in 1:nsig){ set.seed(12545) x = rnorm(N,m0,sqrt(C0)) w = rep(1/N,N) ws = NULL xs = NULL ess = NULL for (t in 1:n){ x1 = rnorm(N,x,tau) w = dnorm(y[t,i],x1,sig[i]) w = w/sum(w) ESS = 1/sum(w^2) x = sample(x1,size=N,replace=T,prob=w) xs = rbind(xs,x) ws = rbind(ws,w/sum(w)) ess = c(ess,ESS) } wss[i,,] = ws xss[i,,] = xs Ess[i,] = ess } mx = matrix(0,n,nsig) Lx = matrix(0,n,nsig) Ux = matrix(0,n,nsig) for (i in 1:nsig){ mx[,i] = apply(xss[i,,],1,median) Lx[,i] = apply(xss[i,,],1,q025) Ux[,i] = apply(xss[i,,],1,q975) } pdf(file=paste("sisr",N,".pdf",sep=""),width=15,height=10) ind = c(1,seq(10,n,by=10)) par(mfrow=c(2,nsig)) for (i in 1:nsig){ L1 = min(L[,i],Lx[,i]) U1 = max(U[,i],Ux[,i]) plot(m[,i],xlab="time",ylab="",ylim=c(L1,U1),lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() lines(L[,i],lwd=2) lines(U[,i],lwd=2) lines(Lx[,i],lwd=2,col=2) lines(mx[,i],lwd=2,col=2) lines(Ux[,i],lwd=2,col=2) } for (i in 1:nsig){ plot(Ess[i,],xlab="time",ylab="",lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() } dev.off() dm = abs((mx-m)) dl = abs((Lx-L)) du = abs((Ux-U)) Lower = quantile(c(dm,dl,du),0.005) Upper = quantile(c(dm,dl,du),0.99) pdf(file=paste("sisr",N,"-abserror.pdf",sep=""),width=15,height=10) par(mfrow=c(1,3)) for (i in 1:nsig) boxplot(dl[,i],dm[,i],du[,i],outline=FALSE,ylim=c(Lower,Upper),names=c("2.5%","50%","97.5%")) dev.off() dl.sisr = dl dm.sisr = dm du.sisr = du # Auxiliary particle filter (APF) # ------------------------------- N = 1000 wss = array(0,c(nsig,n,N)) xss = array(0,c(nsig,n,N)) Ess = matrix(0,nsig,n) for (i in 1:nsig){ set.seed(12545) x = rnorm(N,m0,sqrt(C0)) w = rep(1/N,N) ws = NULL xs = NULL ess = NULL for (t in 1:n){ w0 = dnorm(y[t,i],x,sig[i]) k = sample(1:N,size=N,replace=TRUE,prob=w0) x1 = rnorm(N,x[k],tau) lw = dnorm(y[t,i],x1,sig[i],log=TRUE)-dnorm(y[t,i],x[k],sig[i],log=TRUE) w = exp(lw-max(lw)) w = w/sum(w) ESS = 1/sum(w^2) x = sample(x1,size=N,replace=T,prob=w) xs = rbind(xs,x) ws = rbind(ws,w/sum(w)) ess = c(ess,ESS) } wss[i,,] = ws xss[i,,] = xs Ess[i,] = ess } mx = matrix(0,n,nsig) Lx = matrix(0,n,nsig) Ux = matrix(0,n,nsig) for (i in 1:nsig){ mx[,i] = apply(xss[i,,],1,median) Lx[,i] = apply(xss[i,,],1,q025) Ux[,i] = apply(xss[i,,],1,q975) } pdf(file=paste("apf",N,".pdf",sep=""),width=15,height=10) ind = c(1,seq(10,n,by=10)) par(mfrow=c(2,nsig)) for (i in 1:nsig){ L1 = min(L[,i],Lx[,i]) U1 = max(U[,i],Ux[,i]) plot(m[,i],xlab="time",ylab="",ylim=c(L1,U1),lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() lines(L[,i],lwd=2) lines(U[,i],lwd=2) lines(Lx[,i],lwd=2,col=2) lines(mx[,i],lwd=2,col=2) lines(Ux[,i],lwd=2,col=2) } for (i in 1:nsig){ plot(Ess[i,],xlab="time",ylab="",lwd=2,axes=F,type="l") axis(2);axis(1,at=ind,label=ind);box() } dev.off() dm = abs((mx-m)) dl = abs((Lx-L)) du = abs((Ux-U)) Lower = quantile(c(dm,dl,du),0.005) Upper = quantile(c(dm,dl,du),0.99) pdf(file=paste("apf",N,"-abserror.pdf",sep=""),width=15,height=10) par(mfrow=c(1,3)) for (i in 1:nsig) boxplot(dl[,i],dm[,i],du[,i],outline=FALSE,ylim=c(Lower,Upper),names=c("2.5%","50%","97.5%")) dev.off() dl.apf = dl dm.apf = dm du.apf = du # Comparison # ---------- Lower = quantile(c(dl.sisr,dm.sisr,du.sisr,dl.apf,dm.apf,du.apf),0.005) Upper = quantile(c(dl.sisr,dm.sisr,du.sisr,dl.apf,dm.apf,du.apf),0.8) pdf(file=paste("sisr-apf",N,"-abserror.pdf",sep=""),width=15,height=10) par(mfrow=c(1,3)) for (i in 1:nsig) boxplot(dl.sisr[,i],dl.apf[,i],dm.sisr[,i],dm.apf[,i],du.sisr[,i],du.apf[,i], outline=FALSE,ylim=c(Lower,Upper), names=c("SISR2.5","APF2.5","SISR50","APF50","SISR97.5","APF97.5")) dev.off()