pmhtutorial (version 1.5)

particleFilter: Fully-adapted particle filter for state estimate in a linear Gaussian state space model

Description

Estimates the filtered state and the log-likelihood for a linear Gaussian state space model of the form \( x_{t} = \phi x_{t-1} + \sigma_v v_t \) and \( y_t = x_t + \sigma_e e_t \), where \(v_t\) and \(e_t\) denote independent standard Gaussian random variables, i.e.\(N(0,1)\).

Usage

particleFilter(y, theta, noParticles, initialState)

Arguments

y

Observations from the model for \(t=1,...,T\).

theta

The parameters \(\theta=\{\phi,\sigma_v,\sigma_e\}\) of the LGSS model. The parameter \(\phi\) scales the current state in the state dynamics. The standard deviations of the state process noise and the observation process noise are denoted \(\sigma_v\) and \(\sigma_e\), respectively.

noParticles

The number of particles to use in the filter.

initialState

The initial state.

Value

The function returns a list with the elements:

• xHatFiltered: The estimate of the filtered state at time \(t=1,...,T\).

• logLikelihood: The estimate of the log-likelihood.

• particles: The particle system at each time point.

• weights: The particle weights at each time point.

References

Dahlin, J. & Schon, T. B. "Getting Started with Particle Metropolis-Hastings for Inference in Nonlinear Dynamical Models." Journal of Statistical Software, Code Snippets, 88(2): 1--41, 2019.

Examples

# NOT RUN {
# Generates 500 observations from a linear state space model with
# (phi, sigma_e, sigma_v) = (0.5, 1.0, 0.1) and zero initial state.
theta <- c(0.5, 1.0, 0.1)
d <- generateData(theta, noObservations=500, initialState=0.0) 

# Estimate the filtered state using a Particle filter
pfOutput <- particleFilter(d$y, theta, noParticles = 50, 
  initialState=0.0)

# Plot the estimate and the true state
par(mfrow=c(3, 1))
plot(d$x[1:500], type="l", xlab="time", ylab="true state", bty="n", 
  col="#1B9E77")
plot(pfOutput$xHatFiltered, type="l", xlab="time", 
  ylab="paticle filter estimate", bty="n", col="#D95F02")
plot(d$x[1:500]-pfOutput$xHatFiltered, type="l", xlab="time", 
  ylab="difference", bty="n", col="#7570B3")
# }