A package to perform particle filtering (as well as likelihood estimation and smoothing) using the Feynman-Kac formalism.

Filtering and likelihood estimation are illustrated on two stochastic diffusion equation models:

• The Cox-Ingersoll-Ross (CIR) model
• The K dimensional continuous Wright Fisher model (continuous time, infinite population, see Jenkins & Spanò (2017) for instance)

Particle smoothing for the Wright-Fisher model is not implemented for lack of a tractable form of the transition density.

Outputs:

• Marginal likelihood
• Samples from the filtering distribution
• Samples from the marginal smoothing distribution

Implemented:

• Bootstrap particle filter with adaptive resampling.
• Forward Filtering Backward Sampling (FFBS) algorithm

Potentially useful functions to be found in the package:

• Evaluation of the transition density for the Cox-Ingersoll-Ross process (based on the representation with the Bessel function)
• Random trajectory generation from the Cox-Ingersoll-Ross process (based on the Gamma Poisson expansion of the transition density)

The Feynman-Kac formalism allows to formulate different types of particle filters using the same abstract elements. The input of a generic particle filter are:

• A Feynman-Kac model $M_t, G_t$, where:
  • $G_t$ is a potential function which can be evaluated for all values of $t$
  • It is possible to simulate from $M_0(dx_0)$ and $M_t(x_{t-1}, dx_t)$
• The number of particles $N$
• The choice of an unbiased resampling scheme (e.g. multinomial), i.e. an algorithm to draw variables $A_t^{1:N} \sim RS\left(W_{t-1}^{1:N}\right)$ in $1:N$ where $RS$ is a distribution such that: $\mathbb{E}\left[\sum_{m=1}^N \delta\left(A_t^m=n\right)\right] = W_{t-1}^n$.

For adaptive resampling, one needs in addition:

• a scalar $ESS_{\min} \ge 0$

Using this formalism, the bootstrap filter is expressed as:

• $G_0(x_0) = f_0\bigl(y_0 \mid x_0\bigr)$, where $f$ is the emission density.
• $G_t(x_{t-1}, x_t) = f_t\bigl(y_t \mid x_t\bigr)$ for $t\ge 1$.
• $M_0(dx_0) = P_0(dx_0)$, the prior on the hidden state.
• $M_t(x_{t-1}, dx_t) = P_t\bigl(x_{t-1}, dx_t\bigr)$, given by the transition kernel.

The generic framework of Feynman-Kac models allows to implement particle filtering in an abstract way, independently of the specific nature of the state space or the transition kernel. The potential functions G_t are only required to be evaluable on any state, and the transition kernels M_t to accept a state as input and return another state as output.

Input:

• Mt: sequence of transition simulators (Mt1 ~ M0, Mtt ~ M_t(x,·) for t>=2)
• Gt: sequence of potential functions (Gt1, Gtt) that return non-negative weights
• N: number of particles
• RS(W): resampling routine that returns N ancestor indices given normalized linear weights W
• ESS_min_frac: fraction of N below which to resample (optional)

Output:

• particles[t][i]: particle i at time t
• W[t][i]: normalized (linear) weights at time t
• Z: estimate of marginal likelihood (linear scale)

Procedure:

1. Helpers:

   • normalize(w): W = w / sum(w)
   • ESS(W): 1 / sum(W^2)

2. Initialization (t = 1)

   • For i in 1..N: particles[1][i] = Mt1 # sample from prior M0
   • For i in 1..N: w[i] = Gt1 # non-negative potentials
   • Z = mean(w) # incremental marginal likelihood
   • W[1] = normalize(w)

3. For t = 2..T:

   • Wprev = W[t-1]
   • If ESS(Wprev) < ESS_min_frac * N:
     • ancestors = RS(Wprev) # length N
     • For i in 1..N: particles[t-1][i] = particles[t-1][ancestors[i]]
     • Wprev = [1/N] * N # uniform after resampling
   • For i in 1..N: particles[t][i] = Mtt # propagate
   • For i in 1..N: w[i] = Gt[t](particles[t-1][i], particles[t][i]) # incremental weights
   • W[t] = normalize(w)

4. Return (particles, W, Z)

Press ] in the Julia interpreter to enter the Pkg mode and input:

pkg> add https://github.com/konkam/FeynmanKacParticleFilters.jl

The transition density of the 1-D CIR process is available as:

$$P_t(x, dx') = \sum_{k \ge 0}\text{Poisson}\left(k; \frac{\gamma}{\sigma^2}\frac{1}{e^{2\gamma t}-1}x\right)\text{Ga}\left(x;k+\delta/2, \frac{\gamma}{\sigma^2}\frac{e^{2\gamma t}}{e^{2\gamma t}-1}\right)$$

from which it easy to simulate. Moreover, we consider a Poisson distribution as the emission density:

$$f_t(y_t|x_t) = \frac{x_t^{y_t}}{y_t!}e^{-x_t}$$

We start by simulating some data (a function to simulate from the transition density is available in the package):

using FeynmanKacParticleFilters, Distributions, Random

Random.seed!(0)

Δt = 0.1
δ = 3.
γ = 2.5
σ = 4.
Nobs = 2
Nsteps = 4
λ = 1.
Nparts = 10
α = δ/2
β = γ/σ^2

time_grid = [k*Δt for k in 0:(Nsteps-1)]
times = [k*Δt for k in 0:(Nsteps-1)]
X = FeynmanKacParticleFilters.generate_CIR_trajectory(time_grid, 3, δ*1.2, γ/1.2, σ*0.7)
Y = map(λ -> rand(Poisson(λ), Nobs), X);
data = zip(times, Y) |> Dict

4-element Array{Float64,1}:
 0.0
 0.1
 0.2
 0.30000000000000004

Now we define the (log)potential function Gt, a simulator from the transition kernel for the Cox-Ingersoll-Ross model (we use convenience functions to create all potentials and kernels) and a resampling scheme (here multinomial):

Mt = FeynmanKacParticleFilters.create_transition_kernels_CIR(data, δ, γ, σ)
logGt = FeynmanKacParticleFilters.create_log_potential_functions_CIR(data)
RS(W) = rand(Categorical(W), length(W))

Running the boostrap filter algorithm is performed as follows:

pf = FeynmanKacParticleFilters.generic_particle_filtering_adaptive_resampling_logweights(Mt, logGt, Nparts, RS)

To sample nsamples values from the i-th filtering distributions, do:

n_samples = 100
i = 4
FeynmanKacParticleFilters.sample_from_filtering_distributions_logweights(pf, n_samples, i)
100-element Array{Float64,1}:
  5.371960182098351
  5.371960182098351
  3.3924167451813956
  3.3924167451813956
  3.3924167451813956
  ⋮

To perform a simple particle smoothing on the CIR process using the FFBS algorithm, we additionally need a function which evaluates the transition density of the CIR process.

transition_logdensity_CIR(Xtp1, Xt, Δtp1) = FeynmanKacParticleFilters.CIR_transition_logdensity(Xtp1, Xt, Δtp1, δ, γ, σ)

With the transition density, we can obtain the FFBS filter:

ps = FeynmanKacParticleFilters.generic_FFBS_algorithm_logweights(Mt, logGt, Nparts, Nparts, RS, transition_logdensity_CIR)

To sample nsamples values from the i-th smoothing distribution, do:

n_samples = 100
i = 4
FeynmanKacParticleFilters.sample_from_smoothing_distributions_logweights(ps, n_samples, i)
100-element Array{Float64,1}:
 7.134633585387236
 2.513540876531395
 5.0555536713845814
 7.983322471825221
 4.651221100411266
 ⋮

References:

• Briers, M., Doucet, A. and Maskell, S. Smoothing algorithms for state–space models. Annals of the Institute of Statistical Mathematics 62.1 (2010): 61.

• Chopin, N. & Papaspiliopoulos, O. A concise introduction to Sequential Monte Carlo, to appear.

• Del Moral, P. (2004). Feynman-Kac formulae. Genealogical and interacting particle systems with applications. Probability and its Applications. Springer Verlag, New York.

• Jenkins, P. A., & Spanò, D. (2017). Exact simulation of the Wright--Fisher diffusion. The Annals of Applied Probability, 27(3), 1478–1509.