rm(list=ls()) ffbs = function(y,sig2,tau2,a1,R1,ndraw){ 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) # time t=1 a[1] = a1 R[1] = R1 f = a[1] Q = R[1] + sig2 A = R[1]/Q m[1] = a[1] + A*(y[1]-f) C[1] = R[1] - Q*A**2 # forward filtering for (t in 2:n){ a[t] = m[t-1] R[t] = C[t-1] + tau2 f = a[t] Q = R[t] + sig2 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 } # backward sampling if (ndraw==1){ theta = rep(0,n) theta[n] = rnorm(1,m[n],sqrt(C[n])) for (t in (n-1):1) theta[t] = rnorm(1,m[t]+B[t]*(theta[t+1]-a[t+1]),sqrt(H[t])) }else{ theta = matrix(0,ndraw,n) theta[,n] = rnorm(ndraw,m[n],sqrt(C[n])) for (t in (n-1):1) theta[,t] = rnorm(ndraw,m[t]+B[t]*(theta[,t+1]-a[t+1]),sqrt(H[t])) } return(theta) } sim.dlm = function(n,sig,tau,beta0){ e = rnorm(n,0,sig) u = rnorm(n,0,tau) beta = rep(0,n) y = rep(0,n) beta[1] = beta0 + u[1] for (t in 2:n) beta[t] = beta[t-1] + u[t] for (t in 1:n) y[t] = beta[t] + e[t] return(list(y=y,beta=beta)) } # Computing running times set.seed(154697) sig = 1.0 tau = 0.5 beta0 = 0.0 sig2 = sig^2 tau2 = tau^2 b0 = 0 B0 = 10000 a1 = b0 R1 = B0 + tau2 M0 = 1000 M = 1000 L = 1 niter = M0+L*M ns = c(seq(100,1000,by=100),2000,3000,4000,5000) ii = 0 blrtime = rep(0,length(ns)) ffbstime = rep(0,length(ns)) for (n in ns){ print(n) sim = sim.dlm(n,sig,tau,beta0) y = sim$y blr.time = system.time({ one = matrix(1,n,1) A = diag(1:n) for (i in 1:(n-1)){ A[i,(i+1):n]=i A[(i+1):n,i]=i } iA = solve(A) iAone = iA%*%one V1 = solve(iA/tau2+diag(1/sig2,n)) m1 = V1%*%(beta0*iAone/tau2+y/sig2) }) ffbs.time = system.time({ betas = ffbs(y,sig2,tau2,a1,R1,niter) }) ii = ii + 1 blrtime[ii] = blr.time[3] ffbstime[ii] = ffbs.time[3] } pdf(file="runningtime.pdf",width=8,height=6) cbind(ns,blrtime,ffbstime,blrtime/ffbstime) plot(ns,log(blrtime),xlab="Number of observations",ylab="Running time (LOG)",type="b",lwd=2,axes=FALSE,pch=16) axis(2);box();axis(1,at=ns) lines(ns,log(ffbstime),col=2,type="b",pch=16) legend(100,4,legend=c("BLR","FFBS"),col=1:2,lwd=2,bty="n") dev.off()