Summary

The PF library provides class and function templates that offer fast implementations for a variety of particles filtering algorithms. Each of these algorithms are useful for a wide range of time series models.

In this library, each available particle filtering algorithm is provided as an abstract base class template. Once the data analyst has a specific state-space or hidden Markov model in mind, she will pick which type(s) of particle filtering algorithm(s) to associate with that model by including the appropriate header file. For each model-particle filter pair, she will write a class template for his model that inherits from the particle filter’s base class template.

The details of each algorithm are abstracted away, and each algorithm’s class template’s required functions are pure virtual methods, meaning that the data analyst will not be able to omit any function that is required by the algorithm.

This is by no means the first C++ library to be offered that provides particle filter implementations. Other options include LibBi (Murray 2013), Biips (Todeschini et al. 2014), SMCTC (Johansen 2009) and Nimble (Valpine et al. 2017). The goals of these software packages are different, though–users of these libraries write their models in a scripting language, and that model file gets parsed into C++ code. This library is designed for users that prefer to work in C++ directly.

Statement of Need

It takes time and effort to implement different particle filters well, and this is true for two reasons. First, the mathematical notation used to describe them can be complex. The second reason is that, even if they are correctly implemented, they can be quite slow, limiting the number of tasks that they would be feasible for. This library attempts to provide speed and abstraction to mitigate these two difficulties.

Additionally, this software is designed in an object-oriented manner. This allows for individual particle filters to be used, and it facilitates the implementation of algorithms that update many particle filters, possibly in parallel. Some examples of opportunities include particle filters with parallelized resampling schemes (Bolic, Djuric, and Sangjin Hong 2005, 1309.2918), particle Markov chain Monte Carlo algorithms (Andrieu, Doucet, and Holenstein 2010), importance sampling “squared” (Tran et al. 2013), and the particle swarm algorithm (Brown 2020).

Finally, this library is “header-only.” As a result, building your C++ project is as simple as possible. The only required steps are #include-ing relevant headers, and pointing the compiler at the include/pf/ directory. This directory stores all necessary code, although there there are unit tests and examples provided as well.

Example

A fully-worked example is provided along with this software available at https://github.com/tbrown122387/pf. This example considers modeling a financial time series with a simple stochastic volatility model (Taylor 1982) with three parameters: \(\beta, \phi\) and \(\sigma\). For this model, the observable rate of return \(y_t\) is normally distributed after conditioning on the contemporaneous state random variable \(x_t\). The mean parameter of this normal distribution will remain fixed at \(0\). However, the scale of this distribution will vary with the evolving \(x_t\). When \(x_t\) is relatively high, the returns will have a high conditional variance and be “volatile.” When \(x_t\) is low, the returns will be much less volatile.

The observation equation is \[\begin{eqnarray} y_t = \beta e^{x_t/2} z_t \end{eqnarray}\] where \(\{z_t\}_{t=0}^T\) are independent and identically distributed (iid) normal random variates.

The state evolves randomly through time as an autoregressive process of order one: \[\begin{eqnarray} x_t = \phi x_{t-1} + \sigma z'_t. \end{eqnarray}\] The collection \(\{z'_t\}_{t=1}^T\) are also assumed to be iid normal random variates. At time \(1\), we assume the first state follows a mean zero normal distribution with \(\sigma^2/(1-\phi^2)\). For simplicity, all of our proposal distributions are chosen to be the same as the state transitions.

The file examples/svol_sisr.h provides a fully-worked example of writing a class template called svol_sisr for this model-algorithm pair. Any model-algorithm pair will make use of a resampler type (found in include/pf/resamplers.h), sampler functions (found in include/pf/rv_samp.h), and density evaluator functions (found in include/pf/rv_eval.h). Because you are writing a class template instead of a class, the decision of what to pass in as template parameters will be pushed back to the instantiation site. These template parameters are often changed quite frequently, so this design allows them to be changed only in one location of the project. Instantiating an object after the class template has been written is much easier than writing the class template itself. Just provide the template parameters and the constructor parameters in the correct order. For example,

using Mod = svol_sisr<5000,1,1,mn_resampler<5000,1,double>,double>;
Mod sisrsvol(.91,.5,1.0);

instantiates the object called sisrsvol, which performs the SISR algorithm for the stochastic volatility model using multinomial sampling on \(5,000\) particles. It sets \(\phi = .95\), \(\beta = .5\) and \(\sigma = 1\).

As (possibly real-time, streaming) data becomes available, updating the model is accomplished with the filter method. The call at each time point would look something like

sisrsvol.filter(yt);

References

Andrieu, Christophe, Arnaud Doucet, and Roman Holenstein. 2010. “Particle Markov Chain Monte Carlo Methods.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 (3): 269–342. https://doi.org/10.1111/j.1467-9868.2009.00736.x.

Bolic, M., P. M. Djuric, and Sangjin Hong. 2005. “Resampling Algorithms and Architectures for Distributed Particle Filters.” IEEE Transactions on Signal Processing 53 (7): 2442–50. https://doi.org/10.1109/TSP.2005.849185.

Brown, Taylor R. 2020. “Approximating Posterior Predictive Distributions by Averaging Output from Many Particle Filters.”

Johansen, Adam M. 2009. “SMCTC: Sequential Monte Carlo in C++.” Journal of Statistical Software 30 (6): 1–41. http://www.jstatsoft.org/v30/i06.

Murray, Lawrence M. 2013. “Bayesian State-Space Modelling on High-Performance Hardware Using LibBi.”

Taylor, Stephen. 1982. “Financial Returns Modelled by the Product of Two Stochastic Processes, a Study of Daily Sugar Prices 1961-79” 1 (January).

Todeschini, Adrien, François Caron, Marc Fuentes, Pierrick Legrand, and Pierre Del Moral. 2014. “Biips: Software for Bayesian Inference with Interacting Particle Systems.” arXiv Preprint arXiv:1412.3779.

Tran, Minh-Ngoc, Marcel Scharth, Michael Pitt, and Robert Kohn. 2013. “Importance Sampling Squared for Bayesian Inference in Latent Variable Models.” SSRN Electronic Journal, September. https://doi.org/10.2139/ssrn.2386371.

Valpine, Perry de, Daniel Turek, Christopher J. Paciorek, Clifford Anderson-Bergman, Duncan Temple Lang, and Rastislav Bodik. 2017. “Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE.” Journal of Computational and Graphical Statistics 26 (2): 403–13. https://doi.org/10.1080/10618600.2016.1172487.