FeynmanKacParticleFilters
A package to perform particle filtering (as well as likelihood estimation and smoothing) using the FeynmanKac formalism.
Filtering and likelihood estimation are illustrated on two stochastic diffusion equation models:
 The CoxIngersollRoss (CIR) model
 The K dimensional continuous Wright Fisher model (continuous time, infinite population, see Jenkins & Spanò (2017) for instance)
Particle smoothing for the WrightFisher 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:
 Evaluation of the transition density for the CoxIngersollRoss process (based on the representation with the Bessel function)
 Random trajectory generation from the CoxIngersollRoss process (based on the Gamma Poisson expansion of the transition density)
Preliminary notions on the FeynmanKac formalism
The FeynmanKac formalism allows to formulate different types of particle filters using the same abstract elements. The input of a generic particle filter are:
 A FeynmanKac 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(dx0) and M_t(x_t1, dxt)
 The number of particles N
 The choice of an unbiased resampling scheme (e.g. multinomial), i.e. an algorithm to draw variables in 1:N where RS is a distribution such that: .
For adaptive resampling, one needs in addition:
 a scalar
Using this formalism, the boostrap filter is expressed as:
 G_0(x_0) = f_0(y_0x_0), where f is the emission density
 G_t(x_t1, x_t) = f_0(y_tx_t)
 M_0(dx0) = P_0(dx0) the prior on the hidden state
 M_t(x_t1, dxt) = P_t(x_t1, dxt) given by the transition function
How to install the package
Press ]
in the Julia interpreter to enter the Pkg mode and input:
pkg> add https://github.com/konkam/FeynmanKacParticleFilters.jl
How to use the package (Example with the CIR model)
The transition density of the 1D CIR process is available as:
from which it easy to simulate. Moreover, we consider a Poisson distribution as the emission density:
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:(Nsteps1)]
times = [k*Δt for k in 0:(Nsteps1)]
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
4element Array{Float64,1}:
0.0
0.1
0.2
0.30000000000000004
Filtering
Now we define the (log)potential function Gt, a simulator from the transition kernel for the CoxIngersollRoss model 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 ith filtering distributions, do:
n_samples = 100
i = 4
FeynmanKacParticleFilters.sample_from_filtering_distributions_logweights(pf, n_samples, i)
100element Array{Float64,1}:
5.371960182098351
5.371960182098351
3.3924167451813956
3.3924167451813956
3.3924167451813956
⋮
Smoothing
Forward Filtering Backward Sampling (FFBS)
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 ith smoothing distribution, do:
n_samples = 100
i = 4
FeynmanKacParticleFilters.sample_from_smoothing_distributions_logweights(ps, n_samples, i)
100element 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). FeynmanKac 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 WrightFisher diffusion. The Annals of Applied Probability, 27(3), 1478–1509.