Simple nonlinear dynamic model

We will simulated \(n=100\) observations from the following dynamic model: \[ y_t|x_t \sim t_\nu(x_t,1) \ \ \ \mbox{and} \ \ \ x_t|x_{t-1} \sim N(|x_{t-1}|,1), \] with \(x_0=0\) and \(\nu_4\).

set.seed(31415)
n    = 100
x0   = 0
x    = rep(0,n)
x[1] = rnorm(1,abs(x0),1)
for (t in 2:n)
  x[t] = rnorm(1,abs(x[t-1]),1)
y = x+rt(n,df=4)

plot(y,xlab="time",ylab="")
lines(x,col=2)
legend("topleft",legend=c("y(t)","x(t)"),col=1:2,lty=,lwd=3)

One step particle filter (bootstrap filter)

Below we go from \(\{x_0^{(i)}\}_{i=1}^M \sim p(x_0|y^0) \equiv N(0,10)\) to \(\{{\tilde x}_1^{(i)}\}_{i=1}^M \sim p(x_1|y^0)\) to \(\{x_1^{(i)}\}_{i=1}^M \sim p(x_1|y^1)\), for \(M=10,000\) particles.

M  = 10000
x0 = rnorm(M,0,10)
xs = rnorm(M,abs(x0),1)
w  = dt(y[1]-xs,df=4)
k  = sample(1:M,size=M,replace=TRUE,prob=w)
plot(density(x0))
lines(density(xs),col=2)
lines(density(xs[k]),col=4)
legend("topright",legend=c("Posterior at t=0","Prior at t=1","Posterior at t=1"),col=c(1,2,4),lty=1)

Running the bootstrap filter from \(t=1\) to \(t=n\) with \(M=100,000\).

M   = 100000
xss = matrix(0,n,M)
x0  = rnorm(M,0,10)
xs  = rnorm(M,abs(x0),1)
w   = dt(y[1]-xs,df=4)
k   = sample(1:M,size=M,prob=w,replace=TRUE)
xs  = xs[k]
xss[1,] = xs
for (t in 2:n){
  xs = rnorm(M,abs(xs),1)
  w  = dt(y[t]-xs,df=4)
  k  = sample(1:M,size=M,prob=w,replace=TRUE)
  xs = xs[k]
  xss[t,] = xs
}

par(mfrow=c(2,2))
hist(xss[1,],xlab="",prob=TRUE,ylab="Density",main="p(x(1)|y^(1))")
hist(xss[10,],xlab="",prob=TRUE,ylab="Density",main="p(x(10)|y^(10))")
hist(xss[50,],xlab="",prob=TRUE,ylab="Density",main="p(x(50)|y^(50))")
hist(xss[100,],xlab="",prob=TRUE,ylab="Density",main="p(x(100)|y^(100))")

We compute the particle filter approximation to the 5th, 50th and 95th quantiles of \(p(x_t|y^t)\) for all \(t\).

qpf = t(apply(xss,1,quantile,c(0.05,0.5,0.95)))
par(mfrow=c(1,1))
ts.plot(qpf,xlab="Time",ylab="Filtered densities",main="")