plot.hed = function(m,C,text){ xx = seq(m-4*sqrt(C),m+4*sqrt(C),length=200) plot(xx,dnorm(xx,m,sqrt(C)),xlab="",ylab="",main=text,type="l",yaxt="n",xaxt="n", cex.lab=1.5);axis(2,cex.axis=1.25);axis(1,cex.axis=1.25) } ###################################################### # Simulating data following a 1st order DLM (aka local level model) ###################################################### set.seed(12345) n = 50 V = 0.5 W = 0.25 theta0 = 0.0 v = rnorm(n,0,sqrt(V)) w = rnorm(n,0,sqrt(W)) theta = rep(0,n) y = rep(0,n) theta[1] = theta0+w[1] y[1] = theta[1] + v[1] for (t in 2:n){ theta[t] = theta[t-1]+w[t] y[t] = theta[t] + v[t] } par(mfrow=c(1,2)) plot(y,pch=16,ylim=range(y,theta),xlab="time",ylab="") lines(theta,col=2,type="o",pch=16) legend("topleft",legend=c("y(t)","theta(t)"),col=1:2,lty=1,lwd=2) title("SIMULATED WORLD") plot(y,pch=16,ylim=range(y,theta),xlab="time",ylab="") title("REAL WORLD") ######################## # Kalman filter (step by step) ######################## m0 = 0 C0 = 1 m = m0 C = C0 par(mar=c(5.1,5.1,4.1,2.1)) par(mfrow=c(1,3)) t=1 # Prior at time t a = m R = C+W plot.hed(a,R,paste("p(th[",t,"]|y^[",t-1,"])",sep="")) # One-step ahead forecast f = a Q = R + V plot.hed(f,Q,paste("p(y[",t,"]|y^[",t-1,"])",sep="")) points(y[t],0,col=2,pch=16,cex=2) # Posterior at time t A = R/Q m = a+A*(y[t]-f) C = R-A*A*Q plot.hed(m,C,paste("p(theta[",t,"]|y^[",t,"])",sep="")) ########################################################### # Kalman filter and smoother (implementation of classnotes recursions) ########################################################### m0 = 0 C0 = 100 m = m0 C = C0 mf = rep(0,n) Cf = rep(0,n) af = rep(0,n) Rf = rep(0,n) # Kalman filter for (t in 1:n){ # Prior at time t a = m R = C+W af[t] = a Rf[t] = R # One-step ahead forecast f = a Q = R + V # Posterior at time t A = R/Q m = a+A*(y[t]-f) C = R-A*A*Q # Storing means and variances mf[t] = m Cf[t] = C } # Kalman smoother recursions mb=rep(0,n) Cb=rep(0,n) mb[n]=mf[n] Cb[n]=Cf[n] for (t in (n-1):1){ mb[t] = mf[t] + Cf[t]/Rf[t+1]*(mb[t+1]-af[t+1]) Cb[t] = Cf[t] - Cf[t]/Rf[t+1]*(Rf[t+1]-Cb[t+1])*Cf[t]/Rf[t+1] } # Plotting the results Lf = mf-2*sqrt(Cf) Uf = mf+2*sqrt(Cf) Lb = mb-2*sqrt(Cb) Ub = mb+2*sqrt(Cb) par(mfrow=c(1,1)) plot(mf,xlab="time",ylab="",ylim=range(Lf,Uf,Lb,Ub),type="o",pch=16) lines(Lf) lines(Uf) lines(mb,col=2,type="o",pch=16) lines(Lb,col=2) lines(Ub,col=2) legend("topleft",legend=c("Filtered","Smoothed"),col=1:2,lty=1,lwd=2) ############################################# # Using the R package dlm (Petris and Petrone) ############################################# install.packages("dlm");library(dlm) install.packages("numDeriv");library(numDeriv) install.packages("tseries");library(tseries) install.packages("forecast");library(forecast) dlm0 = function(parm){ return(dlm(FF=1,V=parm[1],GG=1,W=parm[2],m0=0,C0=100)) } dlm0fit = dlm0(c(V,W)) filtered = dlmFilter(y,dlm0fit) smoothed = dlmSmooth(filtered) mf1 = filtered$m[2:(n+1)] Cf1 = filtered$D.C[2:(n+1)]^2 Lf1 = mf1-2*sqrt(Cf1) Uf1 = mf1+2*sqrt(Cf1) mb1 = smoothed$s[2:(n+1)] Cb1 = smoothed$D.S[2:(n+1)]^2 Lb1 = mb1-2*sqrt(Cb1) Ub1 = mb1+2*sqrt(Cb1) par(mfrow=c(1,2)) plot(mf,xlab="time",ylab="",ylim=range(Lf,Uf,Lf1,Uf1),type="l") lines(Lf) lines(Uf) lines(mf1,col=2) lines(Lf1,col=2) lines(Uf1,col=2) title("Filtered") plot(mb,xlab="time",ylab="",ylim=range(Lb,Ub,Lb1,Ub1),type="l") lines(Lb) lines(Ub) lines(mb1,col=2) lines(Lb1,col=2) lines(Ub1,col=2) title("Smoothed") # Obtained final smoothed distribution (after estimating V and W) # ----------------------------------------------------------------------------- fit = dlmMLE(y,parm=c(0.5,0.5),build=dlm0,hessian=T) Vhat=fit$par[1] What=fit$par[2] dlm0fit = dlm0(c(Vhat,What)) filtered = dlmFilter(y,dlm0fit) smoothed = dlmSmooth(filtered) mb = smoothed$s[2:(n+1)] Cb = smoothed$D.S[2:(n+1)]^2 Lb = mb-2*sqrt(Cb) Ub = mb+2*sqrt(Cb) par(mfrow=c(1,1)) plot(mb,xlab="time",ylab=expression(theta[t]),ylim=range(Lb,Ub),type="l") lines(Lb,lty=2) lines(Ub,lty=2) title(paste("Smoothed distributions\nV.hat=",round(Vhat,3)," - W.hat=",round(What,3),sep="")) # Producing h-steps ahead forecasting # ------------------------------------------- h=10 fore = dlmForecast(filtered,nAhead=h,method="plain") yhat = as.numeric(fore$f) Vyhat=as.numeric(fore$Q) Uyhat = yhat+2*sqrt(Vyhat) Lyhat = yhat-2*sqrt(Vyhat) par(mfrow=c(1,1)) plot(y,ylab=expression(y[t]),main="",xlab="Time",pch=16,xlim=c(1,n+h),ylim=range(y,yhat,Lyhat,Uyhat)) points((n+1):(n+h),yhat,col=2,pch=16) lines((n+1):(n+h),Lyhat,col=2) lines((n+1):(n+h),Uyhat,col=2) title(paste(h,"-steps ahead forecast",sep=""))