1. Introduction

Particle filtering, or sequential Monte Carlo (SMC), methodology involves the simulation over time of an artificial particle system (ξⁱ_t; t ∈ {1, …, T}, i ∈ {1, …, N}). It is particularly suited to numerical approximation of integrals of the form (1) $$Z:={\int }_{{\mathsf{X}}^T}{\mu }_1\left(x_1\right)g_1\left(x_1\right)\prod _{t=2}^T{f}_t(x_{t-1},x_t)g_t\left(x_t\right)d{x}_{1:T},$$(1) where $\mathsf{X}={\mathbb{R}}^d$ for some $d\in \mathbb{N}$, $T\in \mathbb{N}$, x_{1: T} ≔ (x₁, …, x_T), μ₁ is a probability density function on $\mathsf{X}$, each f_t a transition density on $\mathsf{X}$, and each g_t is a bounded, continuous, and nonnegative function. Algorithm 1 describes a particle filter, using which an estimate of (1) can be computed as (2) $$Z^N:=\prod _{t=1}^T\left[\frac{1}{N}\sum _{i=1}^N{g}_t\left({\xi }_t^i\right)\right].$$(2)

A Particle Filter

Particle filters were originally applied to statistical inference for hidden Markov models (HMMs) by Gordon, Salmond, and Smith (1993), and this setting remains an important application. Letting $\mathsf{Y}={\mathbb{R}}^{d^{\text{'}}}$ for some $d^{\text{'}}\in \mathbb{N}$, an HMM is a Markov chain evolving on $\mathsf{X}\times \mathsf{Y}$, ${(X_t,Y_t)}_{t\in \mathbb{N}}$, where ${\left(X_t\right)}_{t\in \mathbb{N}}$ is itself a Markov chain and for t ∈ {1, …, T}, each Y_t is conditionally independent of all other random variables given X_t. In a time-homogeneous HMM, letting $\mathbb{P}$ denote the law of this bivariate Markov chain, we have (3) $$\begin{array}{ccc}& & \mathbb{P}(X_{1:T}\in A,Y_{1:T}\in B):={\int }_{A\times B}\mu \left(x_1\right)g(x_1,y_1)\hfill \\ & & \prod _{t=2}^{T}f(x_{t-1},x_t)g(x_t,y_t)d{x}_{1:T}d{y}_{1:T},\hfill \end{array}$$(3) where $\mu :\mathsf{X}\to {\mathbb{R}}_{+}$ is a probability density function, $f:\mathsf{X}\times \mathsf{X}\to {\mathbb{R}}_{+}$ a transition density, $g:\mathsf{X}\times \mathsf{Y}\to {\mathbb{R}}_{+}$ an observation density and A and B measurable subsets of ${\mathsf{X}}^T$ and ${\mathsf{Y}}^T$, respectively. Statistical inference is often conducted upon the basis of a realization y_{1: T} of Y_{1: T} for some finite T, which we will consider to be fixed throughout the remainder of the article. Letting $\mathbb{E}$ denote expectations w.r.t. $\mathbb{P}$, our main statistical quantity of interest is $L:=\mathbb{E}\left[{\prod }_{t=1}^{T}g(X_t,y_t)\right]$, the marginal likelihood associated with y_{1: T}. In the above, we take ${\mathbb{R}}_{+}$ to be the nonnegative real numbers, and assume throughout that L > 0.

Running Algorithm 1 with (4) $${\mu }_1=\mu ,f_t=f,g_t\left(x\right)=g(x,y_t),$$(4) corresponds exactly to running the bootstrap particle filter (BPF) of Gordon, Salmond, and Smith (1993), and we observe that when (4) holds, the quantity Z defined in (1) is identical to L, so that Z^N defined in (2) is an approximation of L. In applications where L is the primary quantity of interest, there is typically an unknown statistical parameter θ ∈ Θ that governs μ, f, and g, and in this setting the map θ↦L(θ) is the likelihood function. We continue to suppress the dependence on θ from the notation until Section 5.

The accuracy of the approximation Z^N has been studied extensively. For example, the expectation of Z^N, under the law of the particle filter, is exactly Z for any $N\in \mathbb{N}$, and Z^N converges almost surely to Z as N → ∞; these can be seen as consequences of Del Moral (2004, Theorem 7.4.2). For practical values of N, however, the quality of the approximation can vary considerably depending on the model and/or observation sequence. When used to facilitate parameter estimation using, for example, particle Markov chain Monte Carlo (Andrieu, Doucet, and Holenstein 2010), it is desirable that the accuracy of Z^N be robust to small changes in the model and this is not typically the case.

In Section 2 we introduce a family of “twisted HMMs,” parameterized by a sequence of positive functions $\mathit{}\mathbf{\psi }:=({\psi }_1,\dots ,{\psi }_T)$. Running a particle filter associated with any of these twisted HMMs provides unbiased and strongly consistent estimates of L. Some specific definitions of $\mathit{}\mathbf{\psi }$ correspond to well-known modifications of the BPF, and the algorithm itself can be viewed as a generalization of the auxiliary particle filter (APF) of Pitt and Shephard (1999). Of particular interest is a sequence ${\mathit{}\mathbf{\psi }}^{*}$ for which Z^N = L with probability 1. In general, ${\mathit{}\mathbf{\psi }}^{*}$ is not known and the corresponding APF cannot be implemented, so our main focus in Section 3 is approximating the sequence ${\mathit{}\mathbf{\psi }}^{*}$ iteratively, and defining final estimates through use of a simple stopping rule. In the applications of Section 5, we find that the resulting estimates significantly outperform the BPF, and exhibit some robustness to both increases in the dimension of the latent state space $\mathsf{X}$ and changes in the model parameters. There are some restrictions on the class of transition densities and the functions ψ₁, …, ψ_T that can be used in practice, which we discuss.

This work builds upon a number of methodological advances, most notably the twisted particle filter (Whiteley and Lee 2014), the APF (Pitt and Shephard 1999), block sampling (Doucet, Briers, and Sénécal 2006), and look-ahead schemes (Lin et al. 2013). In particular, the sequence ${\mathit{}\mathbf{\psi }}^{*}$ is closely related to the generalized eigenfunctions described in Whiteley and Lee (2014), but in that work the particle filter as opposed to the HMM was twisted to define alternative approximations of L. For simplicity, we have presented the BPF in which multinomial resampling occurs at each timestep. Commonly employed modifications of this algorithm include adaptive resampling (Kong, Liu, and Wong 1994; Liu and Chen 1995) and alternative resampling schemes (see, e.g., Douc, Cappé, and Moulines 2005). Generalization to the time-inhomogeneous HMM setting is fairly straightforward, so we restrict ourselves to the time-homogeneous setting for clarity of exposition.

2. Twisted Models and the ψ-Auxiliary Particle Filter

Given an HMM (μ, f, g) and a sequence of observations y_{1: T}, we introduce a family of alternative twisted models based on a sequence of real-valued, bounded, continuous, and positive functions$\psi :=({\psi }_1,{\psi }_2,\cdots ,{\psi }_T)$. Letting, for an arbitrary transition density f and function ψ, $f(x,\psi ):={\int }_{\mathsf{X}}f(x,x^{\text{'}})\psi \left(x^{\text{'}}\right)d{x}^{\text{'}}$, we define a sequence of normalizing functions$({\tilde{\psi }}_1,{\tilde{\psi }}_2,\cdots ,{\tilde{\psi }}_T)$ on $\mathsf{X}$ by ${\tilde{\psi }}_t\left(x_t\right):=f(x_t,{\psi }_{t+1})$ for t ∈ {1, …, T − 1}, ${\tilde{\psi }}_T\equiv 1$, and a normalizing constant ${\tilde{\psi }}_0:={\int }_{\mathsf{X}}\mu \left(x_1\right){\psi }_1\left(x_1\right)d{x}_1$. We then define the twisted model via the following sequence of twisted initial and transition densities: (5) $$\begin{array}{ccc}& & {\mu }_1^{\psi }\left(x_1\right):=\frac{\mu \left(x_1\right){\psi }_1\left(x_1\right)}{{\tilde{\psi }}_0},\hfill \\ & & f_t^{\psi }(x_{t-1},x_t):=\frac{f(x_{t-1},x_t){\psi }_t\left(x_t\right)}{{\tilde{\psi }}_{t-1}\left(x_{t-1}\right)},t\in \{2,\dots ,T\},\hfill \end{array}$$(5) and the sequence of positive functions (6) $$\begin{array}{ccc}& & g_1^{\psi }\left(x_1\right):=g(x_1,y_1)\frac{{\tilde{\psi }}_1\left(x_1\right)}{{\psi }_1\left(x_1\right)}{\tilde{\psi }}_0,\hfill \\ & & g_t^{\psi }\left(x_t\right):=g(x_t,y_t)\frac{{\tilde{\psi }}_t\left(x_t\right)}{{\psi }_t\left(x_t\right)},t\in \{2,\dots ,T\},\hfill \end{array}$$(6) which play the role of observation densities in the twisted model. Our interest in this family is motivated by the following invariance result. We denote by 1 the sequence of constant functions equal to 1 everywhere.

If $\psi$ is a sequence of bounded, continuous and positive functions, and $$Z_{\mathit{}\mathbf{\psi }}:={\int }_{{\mathsf{X}}^T}{\mu }_1^{\mathit{}\mathbf{\psi }}\left(x_1\right)g_1^{\mathit{}\mathbf{\psi }}\left(x_1\right)\prod _{t=2}^T{f}_t^{\mathit{}\mathbf{\psi }}(x_{t-1},x_t)g_t^{\mathit{}\mathbf{\psi }}\left(x_t\right)d{x}_{1:T},$$then $Z_{\mathit{}\mathbf{\psi }}=L$.

We observe that $$\begin{array}{ccc}& & {\mu }_1^{\psi }\left(x_1\right)g_1^{\psi }\left(x_1\right)\prod _{t=2}^T{f}_t^{\psi }\left(x_{t-1},x_t\right)g_t^{\psi }\left(x_t\right)\hfill \\ & & =\frac{\mu \left(x_1\right){\psi }_1\left(x_1\right)}{{\tilde{\psi }}_0}{g}_1\left(x_1\right)\frac{{\tilde{\psi }}_1\left(x_1\right)}{{\psi }_1\left(x_1\right)}{\tilde{\psi }}_0·\hfill \\ & & \prod _{t=2}^T\frac{f\left(x_{t-1},x_t\right){\psi }_t\left(x_t\right)}{{\tilde{\psi }}_{t-1}\left(x_{t-1}\right)}{g}_t^{\mathbf{1}}\left(x_t\right)\frac{{\tilde{\psi }}_t\left(x_t\right)}{{\psi }_t\left(x_t\right)}\hfill \\ & & =\mu \left(x_1\right)g_1^{\mathbf{1}}\left(x_1\right)\prod _{t=2}^{T}f\left(x_{t-1},x_t\right)g_t^{\mathbf{1}}\left(x_t\right),\hfill \end{array}$$and the result follows.

From a methodological perspective, Proposition 1 makes clear a particular sense in which the L.H.S. of (1) is common to an entire family of μ₁, (f_t)_{t ∈ {2, …, T}} and (g_t)_{t ∈ {1, …, T}}. The BPF associated with the twisted model corresponds to choosing (7) $${\mu }_1={\mu }^{\mathit{}\mathbf{\psi }},f_t=f_t^{\mathit{}\mathbf{\psi }},g_t=g_t^{\mathit{}\mathbf{\psi }},$$(7) in Algorithm 1; to emphasize the dependence on $\mathit{}\mathbf{\psi }$, we provide in Algorithm 2 the corresponding algorithm and we will denote approximations of L by $Z_{\mathit{}\mathbf{\psi }}^N$. We demonstrate below that the BPF associated with the twisted model can also be viewed as an APF associated with the sequence $\mathit{}\mathbf{\psi }$, and so refer to this algorithm as the $\mathit{}\mathbf{\psi }$-APF. Since the class of $\mathit{}\mathbf{\psi }$-APF’s is very large, it is natural to consider whether there is an optimal choice of $\mathit{}\mathbf{\psi }$, in terms of the accuracy of the approximation $Z_{\mathit{}\mathbf{\psi }}^N$: the following proposition describes such a sequence.

$\mathit{}\mathbf{\psi }$-Auxiliary Particle Filter

Let ${\mathit{}\mathbf{\psi }}^{*}:=({\psi }_1^{*},\dots ,{\psi }_T^{*})$, where ψ*_T(x_T) ≔ g(x_T, y_T), and (8) $${\psi }_t^{*}\left(x_t\right):=g\left(x_t,y_t\right)\mathbb{E}\left[\prod _{p=t+1}^{T}g\left(X_p,y_p\right)|\left\{X_t=x_t\right\}\right],x_t\in \mathsf{X},$$(8) for t ∈ {1, …, T − 1}. Then, $Z_{{\mathit{}\mathbf{\psi }}^{*}}^N=L$ with probability 1.

It can be established that $$g(x_t,y_t){\tilde{\psi }}_t^{*}\left(x_t\right)={\psi }_t^{*}\left(x_t\right),t\in \{1,\dots ,T\},x_t\in \mathsf{X},$$and so we obtain from (6) that $g_1^{{\psi }^{*}}\equiv {\tilde{\psi }}_0^{*}$ and $g_t^{{\psi }^{*}}\equiv 1$ for t ∈ {2, …, T}. Hence, $$Z_N^{{\psi }^{*}}=\prod _{t=1}^T\left[\frac{1}{N}\sum _{i=1}^N{g}_t^{{\psi }^{*}}\left({\xi }_t^i\right)\right]={\tilde{\psi }}_0^{*},$$with probability 1. To conclude, we observe that $$\begin{array}{ccc}\hfill {\tilde{\psi }}_0^{*}& =& {\int }_{\mathsf{X}}\mu \left(x_1\right){\psi }_1^{*}\left(x_1\right)d{x}_1\hfill \\ & =& {\int }_{\mathsf{X}}\mu \left(x_1\right)\mathbb{E}\left[\prod _{t=1}^{T}g\left(X_t,y_t\right)|\left\{X_1=x_1\right\}\right]d{x}_1\hfill \\ & =& \mathbb{E}\left[\prod _{t=1}^{T}g\left(X_t,y_t\right)\right]=L.\hfill \end{array}$$

Implementation of Algorithm 2 requires that one can sample according to ${\mu }_1^{\mathit{}\mathbf{\psi }}$ and $f_t^{\mathit{}\mathbf{\psi }}(x,·)$ and compute $g_t^{\psi }$ pointwise. This imposes restrictions on the choice of $\mathit{}\mathbf{\psi }$ in practice, since one must be able to compute both ψ_t and ${\tilde{\psi }}_t$ pointwise. In general models, the sequence ${\mathit{}\mathbf{\psi }}^{*}$ cannot be used for this reason as (8) cannot be computed explicitly. However, since Algorithm 2 is valid for any sequence of positive functions $\mathit{}\mathbf{\psi }$, we can interpret Proposition 2 as motivating the effective design of a particle filter by solving a sequence of function approximation problems.

Alternatives to the BPF have been considered before (see, e.g., the “locally optimal” proposal in Doucet, Godsill, and Andrieu 2000 and the discussion in Del Moral 2004, Section 2.4.2). The family of particle filters we have defined using $\mathit{}\mathbf{\psi }$ are unusual, however, in that $g_t^{\mathit{}\mathbf{\psi }}$ is a function only of x_t rather than (x_{t − 1}, x_t); other approaches in which the particles are sampled according to a transition density that is not f typically require this extension of the domain of these functions. This is again a consequence of the fact that the $\mathit{}\mathbf{\psi }$-APF can be viewed as a BPF for a twisted model. This feature is shared by the fully adapted APF of Pitt and Shephard (1999), when recast as a standard particle filter for an alternative model as in Johansen and Doucet (2008), and which is obtained as a special case of Algorithm 2 when ψ_t( · ) ≡ g( ·, y_t) for each t ∈ {1, …, T}. We view the approach here as generalizing that algorithm for this reason.

It is possible to recover other existing methodological approaches as BPFs for twisted models. In particular, when each element of $\mathit{}\mathbf{\psi }$ is a constant function, we recover the standard BPF of Gordon, Salmond, and Smith (1993). Setting ψ_t(x_t) = g(x_t, y_t) gives rise to the fully adapted APF. By taking, for some $k\in \mathbb{N}$ and each t ∈ {1, …, T}, (9) $$\begin{array}{ccc}\hfill {\psi }_t\left(x_t\right)& =& g\left(x_t,y_t\right)\mathbb{E}\left[\prod _{p=t+1}^{(t+k)\wedge T}g\left(X_p,y_p\right)|\left\{X_t=x_t\right\}\right],\hfill \\ & & x_t\in \mathsf{X},\hfill \end{array}$$(9) $\mathit{}\mathbf{\psi }$ corresponds to a sequence of look-ahead functions (see, e.g., Lin et al. 2013) and one can recover idealized versions of the delayed sample method of Chen, Wang, and Liu (2000) (see also the fixed-lag smoothing approach in Clapp and Godsill 1999), and the block sampling particle filter of Doucet, Briers, and Sénécal (2006). When k ⩾ T − 1, we obtain the sequence ${\mathit{}\mathbf{\psi }}^{*}$. Just as ${\mathit{}\mathbf{\psi }}^{*}$ cannot typically be used in practice, neither can the exact look-ahead strategies obtained by using (9) for some fixed k. In such situations, the proposed look-ahead particle filtering strategies are not $\mathit{}\mathbf{\psi }$-APFs, and their relationship to the ${\mathit{}\mathbf{\psi }}^{*}$-APF is consequently less clear. We note that the offline setting we consider here affords us the freedom to define twisted models using the entire data record y_{1: T}. The APF was originally introduced to incorporate a single additional observation, and could therefore be implemented in an online setting, that is, the algorithm could run while the data record was being produced.

3. Function Approximations and the Iterated APF

4. Approximations of Smoothing Expectations

Thus far, we have focused on approximations of the marginal likelihood, L, associated with a particular model and data record y_{1: T}. Particle filters are also used to approximate so-called smoothing expectations, that is, $\pi \left(\phi \right):=\mathbb{E}[\phi \left(X_{1:T}\right)\mid \{Y_{1:T}=y_{1:T}\}]$ for some $\phi :{\mathsf{X}}^{\mathsf{T}}\to \mathbb{R}$. Such approximations can be motivated by a slight extension of (1), $$\begin{array}{c}\hfill \gamma \left(\phi \right):={\int }_{{\mathsf{X}}^T}\phi \left(x_{1:T}\right){\mu }_1\left(x_1\right)g_1\left(x_1\right)\prod _{t=2}^T{f}_t\left(x_{t-1},x_t\right)g_t\left(x_t\right)d{x}_{1:T},\end{array}$$where ϕ is a real-valued, bounded, continuous function. We can write π(ϕ) = γ(ϕ)/γ(1), where 1 denotes the constant function x↦1. We define below a well-known, unbiased, and strongly consistent estimate γ^N(ϕ) of γ(ϕ), which can be obtained from Algorithm 1. A strongly consistent approximation of π(ϕ) can then be defined as γ^N(ϕ)/γ^N(1).

The definition of γ^N(ϕ) is facilitated by a specific implementation of step 2. of Algorithm 1 in which one samples $$\begin{array}{ccc}\hfill A_{t-1}^i& \sim & \mathrm{Categorical}\left(\frac{g_{t-1}\left({\xi }_{t-1}^1\right)}{{\sum }_{j=1}^N{g}_{t-1}\left({\xi }_{t-1}^j\right)},\dots ,\frac{g_{t-1}\left({\xi }_{t-1}^N\right)}{{\sum }_{j=1}^N{g}_{t-1}\left({\xi }_{t-1}^j\right)}\right),\hfill \\ \hfill {\xi }_t^i& \sim & f_t({\xi }_{t-1}^{A_{t-1}^i},·),\hfill \end{array}$$for each i ∈ {1, …, N} independently. Use of, for example, the Alias algorithm (Walker 1974, 1977) gives the algorithm $\mathcal{O}\left(N\right)$ computational complexity, and the random variables (Aⁱ_t; t ∈ {1, …, T − 1}, i ∈ {1, …, N}) provide ancestral information associated with each particle. By defining recursively for each i ∈ {1, …, N}, Bⁱ_T ≔ i and $B_{t-1}^i:=A_{t-1}^{B_t^i}$ for t = T, …, 2, the {1, …, N}^T-valued random variable Bⁱ_{1: T} encodes the ancestral lineage of ξⁱ_T (Andrieu, Doucet, and Holenstein 2010). It follows from Del Moral (2004, Theorem 7.4.2) that the approximation $$\begin{array}{ccc}\hfill {\gamma }^N\left(\phi \right):& =& \left[\frac{1}{N}\sum _{i=1}^N{g}_T\left({\xi }_T^i\right)\phi ({\xi }_1^{B_1^i},{\xi }_2^{B_2^i},\cdots ,{\xi }_T^{B_T^i})\right]\hfill \\ & & \times \prod _{t=1}^{T-1}\left(\frac{1}{N}\sum _{i=1}^N{g}_t\left({\xi }_t^i\right)\right),\hfill \end{array}$$is unbiased and strongly consistent, and a strongly consistent approximation of π(ϕ) is (14) $$\begin{array}{ccc}\hfill {\pi }^N\left(\phi \right):& =& \frac{{\gamma }^N\left(\phi \right)}{{\gamma }^N\left(1\right)}=\frac{1}{{\sum }_{i=1}^N{g}_T\left({\xi }_T^i\right)}\hfill \\ & & \times \sum _{i=1}^N\phi \left({\xi }_1^{B_1^i},{\xi }_2^{B_2^i},\cdots ,{\xi }_T^{B_T^i}\right)g_T\left({\xi }_T^i\right).\hfill \end{array}$$(14) The ${\mathit{}\mathbf{\psi }}^{*}$-APF is optimal in terms of approximating γ(1) ≡ Z and not π(ϕ) for general ϕ. Asymptotic variance expressions akin to Proposition 3, but for ${\pi }_{\mathit{}\mathbf{\psi }}^N\left(\phi \right)$, can be derived using existing results (see, e.g., Del Moral and Guionnet 1999; Chopin 2004; Künsch 2005; Douc and Moulines 2008) in the same manner. These could be used to investigate the influence of $\mathit{}\mathbf{\psi }$ on the accuracy of ${\pi }_{\mathit{}\mathbf{\psi }}^N\left(\phi \right)$ or the interaction between ϕ and the sequence $\mathit{}\mathbf{\psi }$ which minimizes the asymptotic variance of the estimator of its expectation.

Finally, we observe that when the optimal sequence ${\mathit{}\mathbf{\psi }}^{*}$ is used in an APF in conjunction with an adaptive resampling strategy (see Algorithm 5), the weights are all equal, no resampling occurs and the ξⁱ_t are all iid samples from $\mathbb{P}(X_t\in ·\mid \{Y_{1:T}=y_{1:T}\})$. This at least partially justifies the use of iterated $\mathit{}\mathbf{\psi }$-APFs to approximate ${\mathit{}\mathbf{\psi }}^{*}$: the asymptotic variance ${\sigma }_{\mathit{}\mathbf{\psi }}^2$ in (10) is particularly affected by discrepancies between ${\mathit{}\mathbf{\psi }}^{*}$ and $\mathit{}\mathbf{\psi }$ in regions of relatively high conditional probability given the data record y_{1: T}, which is why we have chosen to use the particles as support points to define approximations of ${\mathit{}\mathbf{\psi }}^{*}$ in Algorithm 3.

5. Applications and Examples

The purpose of this section is to demonstrate that the iAPF can provide substantially better estimates of the marginal likelihood L than the BPF at the same computational cost. This is exemplified by its performance when d is large, recalling that $\mathsf{X}={\mathbb{R}}^d$. When d is large, the BPF typically requires a large number of particles to approximate L accurately. In contrast, the ${\mathit{}\mathbf{\psi }}^{*}$-APF computes L exactly, and we investigate below the extent to which the iAPF is able to provide accurate approximations in this setting. Similarly, when there are unknown statistical parameters θ, we show empirically that the accuracy of iAPF approximations of the likelihood L(θ) are more robust to changes in θ than their BPF counterparts.

Unbiased, nonnegative approximations of likelihoods L(θ) are central to the particle marginal Metropolis–Hastings algorithm (PMMH) of Andrieu, Doucet, and Holenstein (2010), a prominent parameter estimation algorithm for general state space hidden Markov models. An instance of a pseudo-marginal Markov chain Monte Carlo algorithm (Beaumont 2003; Andrieu and Roberts 2009), the computational efficiency of PMMH depends, sometimes dramatically, on the quality of the unbiased approximations of L(θ) (Andrieu and Vihola 2015; Lee and Łatuszyński 2014; Sherlock et al. 2015; Doucet et al. 2015) delivered by a particle filter for a range of θ values. The relative robustness of iAPF approximations of L(θ) to changes in θ, mentioned above, motivates their use over BPF approximations in PMMH.

5.2. Linear Gaussian Model

A linear Gaussian HMM is defined by the following initial, transition, and observation Gaussian densities: $\mu (·)=\mathcal{N}(·;m,\Sigma )$, $f(x,·)=\mathcal{N}(·;Ax,B)$, and $g(x,·)=\mathcal{N}(·;Cx,D)$, where $m\in {\mathbb{R}}^d$, $\Sigma ,A,B\in {\mathbb{R}}^{d\times d}$, $C\in {\mathbb{R}}^{d\times d^{\text{'}}}$ and $D\in {\mathbb{R}}^{d^{\text{'}}\times d^{\text{'}}}$. For this model, it is possible to implement the fully adapted APF (FA-APF) and to compute explicitly the marginal likelihood, filtering and smoothing distributions using the Kalman filter, facilitating comparisons. We emphasize that implementation of the FA-APF is possible only for a restricted class of analytically tractable models, while the iAPF methodology is applicable more generally. Nevertheless, the iAPF exhibited better performance than the FA-APF in our examples.

Relative Variance of Approximations of Z When d is Large

We consider a family of linear Gaussian models where m = 0, Σ = B = C = D = I_d and A_ij = α^{|i − j| + 1}, i, j ∈ {1, …, d} for some α ∈ (0, 1). Our first comparison is between the relative errors of the approximations $\widehat{Z}$ of L = Z using the iAPF, the BPF, and the FA-APF. We consider configurations with d ∈ {5, 10, 20, 40, 80} and α = 0.42 and we simulated a sequence of T = 100 observations y_{1: T} for each configuration. We ran 1000 replicates of the three algorithms for each configuration and report box plots of the ratio $\widehat{Z}/Z$ in Figure 1.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 1. Boxplots of $\widehat{Z}/Z$ for different dimensions using 1000 replicates. The crosses indicate the mean of each sample.

Display full size

Figure 1. Boxplots of $\widehat{Z}/Z$ for different dimensions using 1000 replicates. The crosses indicate the mean of each sample.

For^{Figure 2, Figure 3}all the simulations, we ran an iAPF with N₀ = 1000 starting particles, a BPF with N = 10, 000 particles and an FA-APF with N = 5000 particles. The BPF and FA-APF both had slightly larger average computational times than the iAPF with these configurations. The average number of particles for the final iteration of the iAPF was greater than N₀ only in dimensions d = 40 (1033) and d = 80 (1142). For d > 10, it was not possible to obtain reasonable estimates with the BPF in a feasible computational time (similarly for the FA-APF for d > 20). The standard deviation of the samples and the average resampling count across the chosen set of dimensions are reported in Tables 1 and 2.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Table 1. Empirical standard deviation of the quantity $\widehat{Z}/Z$ using 1000 replicates.

Display Table

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Table 2. Average resampling count for the 1000 replicates.

CSVDisplay Table

Fixing the dimension d = 10 and the simulated sequence of observations y_{1: T} with α = 0.42, we now consider the variability of the relative error of the estimates of the marginal likelihood of the observations using the iAPF and the BPF for different values of the parameter α ∈ {0.3, 0.32, …, 0.48, 0.5}. In Figure 2, we report box plots of $\widehat{Z}/Z$ in 1000 replications. For the iAPF, the length of the boxes are significantly less variable across the range of values of α. In this case, we used N = 50, 000 particles for the BPF, giving a computational time at least five times larger than that of the iAPF. This demonstrates that the approximations of the marginal likelihood L(α) provided by the iAPF are relatively insensitive to small changes in α, in contrast to the BPF. Similar simulations, which we do not report, show that the FA-APF for this problem performs slightly worse than the iAPF at double the computational time.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 2. Boxplots of $\widehat{Z}/Z$ for different values of the parameter α using 1000 replicates. The crosses indicate the mean of each sample.

Display full size

Figure 2. Boxplots of $\widehat{Z}/Z$ for different values of the parameter α using 1000 replicates. The crosses indicate the mean of each sample.

Particle Marginal Metropolis–Hastings.

We consider a Linear Gaussian model with m = 0, Σ = B = C = I_d, and D = δI_d with δ = 0.25. We used the lower-triangular matrix $$A=\left(\begin{array}{ccccc}\hfill 0.9& \hfill 0& \hfill 0& 0& 0\\ \hfill 0.3& \hfill 0.7& \hfill 0& 0& 0\\ \hfill 0.1& \hfill 0.2& \hfill 0.6& 0& 0\\ \hfill 0.4& \hfill 0.1& \hfill 0.1& 0.3& 0\\ \hfill 0.1& \hfill 0.2& \hfill 0.5& 0.2& 0\end{array}\right),$$and simulated a sequence of T = 100 observations. Assuming only that A is lower triangular, for identifiability, we performed Bayesian inference for the 15 unknown parameters {A_{i, j}: i, j ∈ {1, …, 5}, j ⩽ i}, assigning each parameter an independent uniform prior on [ − 5, 5]. From the initial point A₁ = I₅, we ran three Markov chains A^BPF_{1: L}, A^iAPF_{1: L}, and A^Kalman_{1: L} of length L = 300, 000 to explore the parameter space, updating one of the 15 parameters components at a time with a Gaussian random walk proposal with variance 0.1. The chains differ in how the acceptance probabilities are computed, and correspond to using unbiased estimates of the marginal likelihood obtain from the BPF, iAPF or the Kalman filter, respectively. In the latter case, this corresponds to running a Metropolis–Hastings (MH) chain by computing the marginal likelihood exactly. We started every run of the iAPF with N₀ = 500 particles. The resulting average number of particles used to compute the final estimate was 500.2. The number of particles N = 20, 000 for the BPF was set to have a greater computational time, in this case A^BPF_{1: L} took 50% more time than A^iAPF_{1: L} to simulate.

In Figure 3, we plot posterior density estimates obtained from the three chains for 3 of the 15 entries of the transition matrix A. The posterior means associated with the entries of the matrix A were fairly close to A itself, the largest discrepancy being around 0.2, and the posterior standard deviations were all around 0.1. A comparison of estimated Markov chain autocorrelations for these same parameters is reported in Figure 4, which indicates little difference between the iAPF-PMMH and Kalman-MH Markov chains, and substantially worse performance for the BPF-PMMH Markov chain. The integrated autocorrelation time of the Markov chains provides a measure of the asymptotic variance of the individual chains’ ergodic averages, and in this regard the iAPF-PMMH and Kalman-MH Markov chains were practically indistinguishable, while the BPF-PMMH performed between 3 and 4 times worse, depending on the parameter. The relative improvement of the iAPF over the BPF does seem empirically to depend on the value of δ. In experiments with larger δ, the improvement was still present but less pronounced than for δ = 0.25. We note that in this example, ${\mathit{}\mathbf{\psi }}^{*}$ is outside the class of possible $\mathit{}\mathbf{\psi }$ sequences that can be obtained using the iAPF: the approximations in $\mathit{}\mathbf{\Psi }$ are functions that are constants plus a multivariate normal density with a diagonal covariance matrix, while the functions in ${\mathit{}\mathbf{\psi }}^{*}$ are multivariate normal densities whose covariance matrices have nonzero, off-diagonal entries.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 3. Linear Gaussian model: density estimates for the specified parameters from the three Markov chains.

Display full size

Figure 3. Linear Gaussian model: density estimates for the specified parameters from the three Markov chains.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 4. Linear Gaussian model: autocorrelation function estimates for the BPF-PMMH (crosses), iAPF-PMMH (solid lines), and Kalman-MH (circles) Markov chains.

Display full size

Figure 4. Linear Gaussian model: autocorrelation function estimates for the BPF-PMMH (crosses), iAPF-PMMH (solid lines), and Kalman-MH (circles) Markov chains.

5.3. Univariate Stochastic Volatility Model

A simple stochastic volatility model is defined by $\mu (·)=\mathcal{N}(·;0,{\sigma }^2/{(1-\alpha )}^2)$, $f(x,·)=\mathcal{N}(·;\alpha x,{\sigma }^2)$ and $g(x,·)=\mathcal{N}(·;0,{\beta }^{2}exp\left(x\right))$, where α ∈ (0, 1), β > 0 and σ² > 0 are statistical parameters (see, e.g., Kim, Shephard, and Chib 1998). To compare the efficiency of the iAPF and the BPF within a PMMH algorithm, we analyzed a sequence of T = 945 observations y_{1: T}, which are mean-corrected daily returns computed from weekday close exchange rates r_{1: T + 1} for the pound/dollar from 1/10/81 to 28/6/85. These data have been previously analyzed using different approaches, for example, in Harvey, Ruiz, and Shephard (1994) and Kim, Shephard, and Chib (1998).

We wish to infer the model parameters θ = (α, σ, β) using a PMMH algorithm and compare the two cases, where the marginal likelihood estimates are obtained using the iAPF and the BPF. We placed independent inverse Gamma prior distributions $\mathcal{IG}(2.5,0.025)$ and $\mathcal{IG}(3,1)$ on σ² and β², respectively, and an independent Beta(20, 1.5) prior distribution on the transition coefficient α. We used $({\alpha }_0,{\sigma }_0,{\beta }_0)=(0.95,\sqrt{0.02},0.5)$ as the starting point of the three chains: X^iAPF_{1: L}, X^BPF_{1: L} and $X_{L^{\text{'}}}^{{\mathrm{BPF}}^{\text{'}}}$. All the chains updated one component at a time with a Gaussian random walk proposal with variances (0.02, 0.05, 0.1) for the parameters (α, σ, β). X^iAPF_{1: L} has a total length of L = 150, 000 and for the estimates of the marginal likelihood that appear in the acceptance probability we use the iAPF with N₀ = 100 starting particles. For X^BPF_{1: L} and $X_{1:L^{\text{'}}}^{{\mathrm{BPF}}^{\text{'}}}$ we use BPFs: X^BPF′_{1: L} is a shorter chain with more particles (L = 150, 000 and N = 1000), while $X_{1:L^{\text{'}}}^{{\mathrm{BPF}}^{\text{'}}}$ is a longer chain with fewer particles (L = 1, 500, 000, N = 100). All chains required similar running time overall to simulate. Figure 5 shows estimated marginal posterior densities for the three parameters using the different chains.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 5. Stochastic volatility model: PMMH density estimates for each parameter from the three chains.

Display full size

Figure 5. Stochastic volatility model: PMMH density estimates for each parameter from the three chains.

In Table 3, we provide the adjusted sample size of the Markov chains associated with each of the parameters, obtained by dividing the length of the chain by the estimated integrated autocorrelation time associated with each parameter. We can see an improvement using the iAPF, although we note that the BPF-PMMH algorithm appears to be fairly robust to the variability of the marginal likelihood estimates in this particular application.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Table 3. Sample size adjusted for autocorrelation for each parameter from the three chains.

CSVDisplay Table

Since particle filters provide approximations of the marginal likelihood in HMMs, the iAPF can also be used in alternative parameter estimation procedures, such as simulated maximum likelihood (Lerman and Manski 1981; Diggle and Gratton 1984). The use of particle filters for approximate maximum likelihood estimation (see, e.g., Kitagawa 1998; Hürzeler and Künsch 2001) has recently been used to fit macroeconomic models (Fernández-Villaverde and Rubio-Ramírez 2007). In Figure 6, we show the variability of the BPF and iAPF estimates of the marginal likelihood at points in a neighborhood of the approximate MLE of (α, σ, β) = (0.984, 0.145, 0.69). The iAPF with N₀ = 100 particles used 100 particles in the final iteration to compute the likelihood in all simulations, and took slightly more time than the BPF with N = 1000 particles, but far less time than the BPF with N = 10, 000 particles. The results indicate that the iAPF estimates are significantly less variable than their BPF counterparts and may therefore be more suitable in simulated maximum likelihood approximations.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 6. Log-likelihood estimates in a neighborhood of the MLE. Boxplots correspond to 100 estimates at each parameter value given by three particle filters, from left to right: BPF (N = 1000), BPF (N = 10, 000), iAPF (N₀ = 100).

Display full size

Figure 6. Log-likelihood estimates in a neighborhood of the MLE. Boxplots correspond to 100 estimates at each parameter value given by three particle filters, from left to right: BPF (N = 1000), BPF (N = 10, 000), iAPF (N₀ = 100).

5.4. Multivariate Stochastic Volatility Model

We consider a version of the multivariate stochastic volatility model defined for $\mathsf{X}={\mathbb{R}}^d$ by $\mu (·)=\mathcal{N}(·;m,U_{★})$, $f(x,·)=\mathcal{N}(·;m+\mathrm{diag}(\varphi \left)\right(x-m),U)$ and $g(x,·)=\mathcal{N}(·;0,exp(\mathrm{diag}\left(x\right)\left)\right)$, where $m,\varphi \in {\mathbb{R}}^d$ and the covariance matrix $U\in {\mathbb{R}}^{d\times d}$ are statistical parameters. The matrix U_⋆ is the stationary covariance matrix associated with (φ, U). This is the basic MSV model in Chib, Omori, and Asai (2009, Sec. 2), with the exception that we consider a nondiagonal transition covariance matrix U and a diagonal observation matrix.

We analyzed two 20-dimensional sequences of observations y_{1: T} and $y_{1:T^{\text{'}}}^{\text{'}}$, where T = 102 and T′ = 90. The sequences correspond to the monthly returns for the exchange rate with respect to the US dollar of a range of 20 different international currencies, in the periods 3/2000–8/2008 (y_{1: T}, pre-crisis) and 9/2008–2/2016 ($y_{1:T^{\text{'}}}^{\text{'}}$, post-crisis), as reported by the Federal Reserve System (available at http://www.federalreserve.gov/releases/h10/hist/). We infer the model parameters θ = (m, φ, U) using the iAPF to obtain marginal likelihood estimates within a PMMH algorithm. A similar study using a different approach and with a set of six currencies can be found in Liu and West (2001).

The aim of this study is to showcase the potential of the iAPF in a scenario where, due to the relatively high dimensionality of the state space, the BPF systematically fails to provide reasonable marginal likelihood estimates in a feasible computational time. To reduce the dimensionality of the parameter space we consider a band diagonal covariance matrix U with nonzero entries on the main, upper, and lower diagonals. We placed independent inverse Gamma prior distributions with mean 0.2 and unit variance on each entry of the diagonal of U, and independent symmetric triangular prior distributions on [ − 1, 1] on the correlation coefficients $\rho \in {\mathbb{R}}^{19}$ corresponding to the upper and lower diagonal entries. We place independent Uniform(0, 1) prior distributions on each component of φ and an improper, constant prior density for m. This results in a 79-dimensional parameter space. As the starting point of the chains, we used φ₀ = 0.95 · 1, diag(U₀) = 0.2 · 1 and for the 19 correlation coefficients we set ρ₀ = 0.25 · 1, where 1 denotes a vector of 1s whose length can be determined by context. Each entry of m₀ corresponds to the logarithm of the standard deviation of the observation sequence of the relative currency.

We ran two Markov chains X_{1: L} and X′_{1: L}, corresponding to the data sequences y_{1: T} and $y_{1:T^{\text{'}}}^{\text{'}}$, both of them updated one component at a time with a Gaussian random walk proposal with standard deviations (0.2 · 1, 0.005 · 1, 0.02 ·  1, 0.02 · 1) for the parameters (m, φ, diag(U), ρ). The total number of updates for each parameter is L = 12, 000 and the iAPF with N₀ = 500 starting particles is used to estimate marginal likelihoods within the PMMH algorithm. In Figure 7, we report the estimated smoothed posterior densities corresponding to the parameters for the Pound Sterling/U.S. Dollar exchange rate series. Most of the posterior densities are different from their respective prior densities, and we also observe qualitative differences between the pre- and post-crisis regimes. For the same parameters, sample sizes adjusted for autocorrelation are reported in Table 4. Considering the high-dimensional state and parameter spaces, these are satisfactory. In the later steps of the PMMH chain, we recorded an average number of iterations for the iAPF of around 5 and an average number of particles in the final $\mathit{}\mathbf{\psi }$-APF of around 502. To compare this with the BPF, we found that with N = 10⁷ particles, at a cost of roughly 200 times the iAPF-PMMH per step, a BPF-PMMH Markov chain accepted only 4 out of 1000 proposed moves, each of which attempted to update a single parameter.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 7. Multivariate stochastic volatility model: density estimates for the parameters related to the Pound Sterling. Pre-crisis chain (solid line), post-crisis chain (dashed line), and prior density (dotted line). The prior densities for (a) and (b) are constant.

Display full size

Figure 7. Multivariate stochastic volatility model: density estimates for the parameters related to the Pound Sterling. Pre-crisis chain (solid line), post-crisis chain (dashed line), and prior density (dotted line). The prior densities for (a) and (b) are constant.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Table 4. Sample size adjusted for autocorrelation.

CSVDisplay Table

The aforementioned qualitative change of regime seems to be evident looking at the difference between the posterior expectations of the parameter m for the post-crisis and the pre-crisis chain, reported in Figure 8. The parameter m can be interpreted as the period average of the mean-reverting latent process of the log-volatilities for the exchange rate series. Positive values of the differences for close to all of the currencies suggest a generally higher volatility during the post-crisis period.

The Iterated Auxiliary Particle Filter

All authors

Pieralberto Guarniero, Adam M. Johansen & Anthony Lee

https://doi.org/10.1080/01621459.2016.1222291

Published online:

18 July 2017

Figure 8. Multivariate stochastic volatility model: differences between post-crisis and pre-crisis posterior expectation of the parameter m for the 20 currencies.

Display full size

Figure 8. Multivariate stochastic volatility model: differences between post-crisis and pre-crisis posterior expectation of the parameter m for the 20 currencies.

6. Discussion

In this article, we have presented the iAPF, an offline algorithm that approximates an idealized particle filter, whose marginal likelihood estimates have zero variance. The main idea is to iteratively approximate a particular sequence of functions, and an empirical study with an implementation using parametric optimization for models with Gaussian transitions showed reasonable performance in some regimes for which the BPF was not able to provide adequate approximations. We applied the iAPF to Bayesian parameter estimation in general state-space HMMs by using it as an ingredient in a PMMH Markov chain. It could also conceivably be used in similar, but inexact, noisy Markov chains; Medina-Aguayo, Lee, and Roberts (2015) showed that control on the quality of the marginal likelihood estimates can provide theoretical guarantees on the behavior of the noisy Markov chain. The performance of the iAPF marginal likelihood estimates also suggests they may be useful in simulated maximum likelihood procedures. In our empirical studies, the number of particles used by the iAPF was orders of magnitude smaller than would be required by the BPF for similar approximation accuracy, which may be relevant for models in which space complexity is an issue.

In the context of likelihood estimation, the perspective brought by viewing the design of particle filters as essentially a function approximation problem has the potential to significantly improve the performance of such methods in a variety of settings. There are, however, a number of alternatives to the parametric optimization approach described in Section 5.1, and it would be of particular future interest to investigate more sophisticated schemes for estimating ${\mathit{}\mathbf{\psi }}^{*}$, that is, specific implementations of Algorithm 3. We have used nonparametric estimates of the sequence ${\mathit{}\mathbf{\psi }}^{*}$ with some success, but the computational cost of the approach was much larger than the parametric approach. Alternatives to the classes $\mathcal{F}$ and Ψ described in Section 3.2 could be obtained using other conjugate families, (see, e.g., Vidoni 1999). We also note that although we restricted the matrix Σ in (15) to be diagonal in our examples, the resulting iAPF marginal likelihood estimators performed fairly well in some situations, where the optimal sequence ${\mathit{}\mathbf{\psi }}^{*}$ contained functions that could not be perfectly approximated using any function in the corresponding class. Finally, the stopping rule in the iAPF, described in Algorithm 4 and which requires multiple independent marginal likelihood estimates, could be replaced with a stopping rule based on the variance estimators proposed in Lee and Whiteley (2015). For simplicity, we have discussed particle filters in which multinomial resampling is used; a variety of other resampling strategies (see Douc, Cappé, and Moulines 2005, for a review) can be used instead.