{'fraction_non_alphanumeric': 0.052728226652675764, 'fraction_numerical': 0.02053456919668882, 'mean_word_length': 4.205750701767696, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Particle Methods for Parameter Estimation in State-Space Models

Nonlinear non-Gaussian state-space models are ubiquitous in statistics, econometrics, information engineering and signal processing. Particle methods, also known as Sequential Monte Carlo (SMC) methods, provide reliable numerical approximations to the associated state inference problems. However, in most applications, the state-space model of interest also depends on unknown static parameters that need to be estimated from the data. In this context, standard particle methods fail and it is necessary to rely on more sophisticated algorithms. The aim of this paper is to present a comprehensive review of particle methods that have been proposed to perform static parameter estimation in state-space models. We discuss the advantages and limitations of these methods and illustrate their performance on simple models.


INTRODUCTION
State-space models, also known as hidden Markov models, are a very popular class of time series models that have found numerous of applications in fields as diverse as statistics, ecology, econometrics, engineering and environmental sciences; see [10], [32], [35], [92]. Formally a statespace model is defined by two stochastic processes {X n } n≥0 and {Y n } n≥0 . The process {X n } n≥0 is a X -valued latent Markov process of initial density µ θ (x) and Markov transition density f θ (x ′ |x), that is whereas the Y-valued observations {Y n } n≥0 satisfy (1.2) Y n | (X 0:n = x 0:n , Y 0:n−1 = y 0:n−1 ) ∼ g θ (y n |x n ),  1 where g θ (y|x) denotes the conditional marginal density, θ ∈ Θ the parameter of the model and z i:j denotes components (z i , z i+1 , ..., z j ) of a sequence {z n }. The spaces X and Y can be Euclidean but what follows applies to more general state spaces as well.
The popularity of state-space models stems from the fact that they are flexible and easily interpretable. Applications of state-space models include stochastic volatility models where X n is the volatility of an asset and Y n its observed log-return [56], biochemical network models where X n corresponds to the population of various biochemical species and Y n are imprecise measurements of the size of a subset of these species [94], neuroscience models where X n is a state vector determining the neuron's stimulus-response function and Y n some spike train data [81]. However, nonlinear non-Gaussian state-space models are also notoriously difficult to fit to data and it is only recently, thanks to the advent of powerful simulation techniques, that it has been possible to fully realize their potential.
To illustrate the complexity of inference in state-space models, consider first the scenario where the parameter θ is known. On-line and off-line inference about the state process {X n } given the observations {Y n } is only feasible analytically for simple models such as the linear Gaussian statespace model. In nonlinear non-Gaussian scenarios, numerous approximation schemes such as the Extended Kalman filter or the Gaussian sum filter [1] have been proposed over the past fifty years to solve these so-called optimal filtering and smoothing problems, but these methods lack rigor and can be unreliable in practice in terms of accuracy, while deterministic integration methods are difficult to implement. Markov chain Monte Carlo (MCMC) methods can obviously be used but they are impractical for on-line inference; and even for off-line inference it can be difficult to build efficient high-dimensional proposal distributions for such algorithms. For nonlinear non-Gaussian state space models particle algorithms have emerged as the most successful. Their widespread popularity is due to the fact that they are easy to implement, suitable for parallel implementation [64] and, more importantly, have been demonstrated in numerous settings to yield more accurate estimates than the standard alternatives; e.g. see [10], [25], [32], [71].
In most practical situations, the model (1.1)-(1.2) depends on an unknown parameter vector θ that needs to be inferred from the data either in an on-line or off-line manner. In fact inferring the parameter θ is often the primary problem of interest; e.g. for biochemical networks, we are not interested in the population of the species per se, but we want to infer some chemical rate constants, which are parameters of the transition prior f θ (x ′ |x). Although it is possible to define an extended state that includes the original state X n and the parameter θ and then apply standard particle methods to perform parameter inference, it was recognized very early on that this naive approach is problematic [58] due to the parameter space not being explored adequately. This has motivated over the past fifteen years the development of many particle methods for the parameter estimation problem, but numerically robust methods have only been proposed recently. The main objective of this paper is to provide a comprehensive overview of this literature. This paper thus differs from recent survey papers on particle methods which all primarily focus on estimating the state sequence X 0:n or discuss a much wider range of topics, e.g. [33], [60], [62], [69]. We will present the main features of each method and comment on their pros and cons. No attempt however is made to discuss the intricacies of the specific implementations. For this we refer the reader to the original references.
We have chosen to broadly classify the methods as follows: Bayesian or Maximum Likelihood (ML) and whether they are implemented off-line or on-line. In the Bayesian approach, the unknown parameter is assigned a prior distribution and the posterior density of this parameter given the observations is to be characterized. In the ML approach, the parameter estimate is the maximizing argument of the likelihood of θ given the data. Both these inference procedures can be carried out offline or on-line. Specifically, in an off-line framework we infer θ using a fixed observation record y 0:T . In contrast, on-line methods update the parameter estimate sequentially as observations {y n } n≥0 become available.
The rest of the paper is organized as follows. In Section 2, we present the main computational challenges associated to parameter inference in state-space models. In Section 3, we review particle methods for filtering when the model does not include any unknown parameters whereas Section 4 is dedicated to smoothing. These filtering and smoothing techniques are at the core of the off-line and on-line ML parameter procedures described in Section 5. In Section 6, we discuss particle methods for off-line and on-line Bayesian parameter inference. The performance of some of these algorithms are illustrated on simple examples in Section 7. Finally, we summarize the main advantages and drawbacks of the methods presented and discuss some open problems in Section 8.

COMPUTATIONAL CHALLENGES ASSOCIATED TO PARAMETER INFERENCE
A key ingredient of ML and Bayesian parameter inference is the likelihood function p θ (y 0:n ) of θ which satisfies where p θ (x 0:n , y 0:n ) denotes the joint density of (X 0:n , Y 0:n ) which is given from equations (1.1)-(1.2) by The likelihood function is also the normalizing constant of the posterior density p θ (x 0:n |y 0:n ) of the latent states X 0:n given data y 0:n (2.3) p θ (x 0:n |y 0:n ) = p θ (x 0:n , y 0:n ) p θ (y 0:n ) .
The main practical issue associated to parameter inference in nonlinear non-Gaussian state-space models is that the likelihood function is intractable. As performing ML parameter inference requires maximizing this intractable function, it means practically that it is necessary to obtain reasonably low variance Monte Carlo estimates of it, or of the associated score vector if this maximisation is carried out using gradient-based methods. Both tasks involve approximating high dimensional integrals, (2.1) and (2.4), whenever n is large. On-line inference requires additionally these integrals to be approximated on the fly, ruling out the applications of standard computational tools such as MCMC.
Bayesian parameter inference is even more challenging, as it requires approximating the posterior density p(θ|y 0:n ) = p θ (y 0:n )p(θ) p θ (y 0:n )p(θ)dθ , (2.5) where p(θ) is the prior density. Here not only p θ (y 0:n ) but also p (y 0:n ) =´p θ (y 0:n )p(θ)dθ are intractable and, once more, these integrals must be approximated on-line if one wants to update the posterior density sequentially. We will show in this review that particle methods are particularly well-suited to these integration tasks.

FILTERING AND PARTICLE APPROXIMATIONS
In this section, the parameter θ is assumed known and we focus on the problem of estimating the latent process {X n } n≥0 sequentially given the observations. An important byproduct of this so-called filtering task from a parameter estimation viewpoint is that it provides us with an on-line scheme to compute {p θ (y 0:n )} n≥0 . As outlined in Section 2, the particle approximation of these likelihood terms is a key ingredient of numerous particle-based parameter inference techniques discussed further on.

Filtering
Filtering denotes usually the task of estimating recursively in time the sequence of marginal posteriors {p θ ( x n | y 0:n )} n≥0 , known as the filtering densities. However we will adopt here a more general definition and will refer to filtering as the task of estimating the sequence of joint posteriors {p θ ( x 0:n | y 0:n )} n≥0 recursively in time but we will still refer to the marginals {p θ ( x n | y 0:n )} n≥0 as the filtering densities.
It is easy to verify from (2.1) and (2.3) that the posterior p θ (x 0:n |y 0:n ) and the likelihood p θ (y 0:n ) satisfy the following fundamental recursions: for n ≥ 1, and (3.2) p θ (y 0:n ) = p θ (y 0:n−1 ) p θ ( y n | y 0:n−1 ) where There are essentially two classes of models for which p θ (x 0:n |y 0:n ) and p θ (y 0:n ) can be computed exactly: the class of linear Gaussian models, for which the above recursions may be implemented using Kalman techniques, and when X is a finite state-space; see for example [10]. For other models these quantities are typically intractable, i.e. the densities in (3.1)-(3.3) cannot be computed exactly.

Algorithm
Particle filtering methods are a set of simulation-based techniques which approximate numerically the recursions (3.1) to (3.3). We focus here on the APF (auxiliary particle filter [82]) for two reasons: first, this is a popular approach, in particular in the context of parameter estimation (see e.g. Section 6.2.3); second, the APF covers as special cases a large class of particle algorithms, such as the bootstrap filter [49] and SISR (Sequential Importance Sampling Resampling [31], [73]). Let (3.4) q θ (x n , y n |x n−1 ) = q θ (x n |y n , x n−1 ) q θ (y n |x n−1 ) where q θ (x n |y n , x n−1 ) is a probability density function which is easy to sample from, and q θ (y n |x n−1 ) is not necessarily required to be a probability density function but just an non-negative function of (x n−1 , y n ) ∈ X × Y one can evaluate. (For n = 0, remove the dependency on x n−1 , i.e. q θ (x 0 , y 0 ) = q θ (x 0 |y 0 )q θ (y 0 ).) The algorithm relies on the following importance weights In order to alleviate the notational burden we omit the dependence of the importance weights on θ; we will do so in the remainder of the paper when no confusion is possible. The auxiliary particle filter can be summarized as follows.

Algorithm 1 Auxiliary Particle Filtering
• At time n = 0 , for all i ∈ {1, ..., N }: • At time n ≥ 1 , for all i ∈ {1, ..., N }: One recovers the SISR algorithm as a special case of Algorithm 1 by taking q θ (y n |x n−1 ) = 1 (or more generally, by taking q θ (y n |x n−1 ) = h θ (y n ), some arbitrary positive function). Further, one recovers the bootstrap filter by taking q θ (x n |y n , x n−1 ) = f θ (x n |x n−1 ). This is an important special case, as some complex models are such that one may sample from f θ (x n |x n−1 ), but not compute the corresponding density; in such a case the bootstrap filter is the only implementable algorithm. For models such that the density f θ (x n |x n−1 ) is tractable, [82] recommend selecting q θ (x n |y n , x n−1 ) = p θ (x n |y n , x n−1 ) and q θ (y n |x n−1 ) = p θ (y n |x n−1 ) when these quantities are tractable, and using approximations of these quantities in scenarios when they are not. The intuition for these recommendations is that this should make the weight function (3.6) nearly constant.
The computational complexity of Algorithm 1 is O (N ) per time step; in particular see e.g. [31, p. 201] for a O (N ) implementation of the resampling step. At time n, the approximations of p θ ( x 0:n | y 0:n ) and p θ ( y n | y 0:n−1 ) presented earlier in (2.3) and (3.3) respectively are given by In practice, one uses (3.7) mostly to obtain approximations of posterior moments N i=1 W i n ϕ(X i 0:n ) ≈ E [ϕ(X 0:n )|y 0:n ] but expressing particle filtering as a method for approximating distributions (rather than moments) turns out to be a more convenient formalization. The likelihood (3.2) is then estimated through The resampling procedure is introduced to replicate particles with high weights and discard particles with low weights. It serves to focus the computational efforts on the "promising" regions of the state-space. We have presented above the simplest resampling scheme. Lower variance resampling schemes have been proposed in [57], [73], as well as more advanced particle algorithms with better overall performance, e.g. the Resample-Move algorithm [47]. For the sake of simplicity, we have also presented a version of the algorithm that operates resampling at every iteration n. It may be more efficient to trigger resampling only when a certain criterion regarding the degeneracy of the weights is met; see [31], [72, pages 35 and 74].

Convergence results
Many sharp convergence results are available for particle methods [25]. A selection of these results that gives useful insights on the difficulties of estimating static parameters with particle methods is presented below.
Under minor regularity assumptions, one can show that for any n ≥ 0, N > 1 and any bounded test function ϕ n : X n+1 → [−1, 1], there exist constants A θ,n,p < ∞ such that for any p ≥ 1 where the expectation is with respect to the law of the particle filter. In addition, for more general classes of functions, we can obtain for any fixed n a Central Limit Theorem (CLT) as N → +∞ [18], [25,Proposition 9.4.2]. Such results are reassuring but weak as they reveal nothing regarding long time behavior. For instance, without further restrictions on the class of functions ϕ n and the state-space model, A θ,n,p typically grows exponentially with n. This is intuitively not surprising as the dimension of the target density p θ ( x 0:n | y 0:n ) is increasing with n. Moreover the successive resampling steps lead to a depletion of the particle population; p θ ( x 0:m | y 0:n ) will eventually be approximated by a single unique particle as n − m increases. This is referred to as the degeneracy problem in the literature [10, Figure 8.4, p. 282]. This is a fundamental weakness of particle methods: given a fixed number of particles N , it is impossible to approximate p θ (x 0:n |y 0:n ) accurately when n is large enough. Fortunately, it is also possible to establish much more positive results. Many state-space models possess the so-called exponential forgetting property [25,Chapter 4]. This property states that for any x 0 , x ′ 0 ∈ X and data y 0:n , there exist constants B θ < ∞ and λ ∈ [0, 1) such that where · TV is the total variation distance; that is the optimal filter forgets exponentially fast its initial condition. This property is typically satisfied when the signal process {X n } n≥0 is a uniformly 6 ergodic Markov chain and the observations {Y n } n≥0 are not too informative [25,Chapter 4] or when {Y n } n≥0 are informative enough that it effectively restricts the hidden state to a bounded region around it [80]. Weaker conditions can be found in [30,97]. When exponential forgetting holds, it is possible to establish much stronger uniform-in-time convergence results for functions ϕ n that depend only on recent states. Specifically, for an integer L > 0 and any bounded test function Ψ L : X L → [−1, 1], there exist constants C θ,L,p < ∞ such that for any p ≥ 1, n ≥ L − 1, This result explains why particle filtering is an effective computational tool in many applications such as tracking, where one is only interested in p θ ( x n−L+1:n | y 0:n ), as the approximation error is uniformly bounded over time. Similar positive results holds for p θ (y 0:n ). This estimate is unbiased for any N ≥ 1 [25, Theorem 7.4.2, page 239] and, under assumption (3.11), the relative variance of the likelihood estimate p θ (y 0:n ), that is the variance of the ratio p θ (y 0:n ) /p θ (y 0:n ), is bounded above by D θ n/N [15], [97]. This is a great improvement over the exponential increase with n that holds for standard importance sampling techniques, see for instance [33]. However, the constants C θ,L,p and D θ are typically exponential in n x , the dimension of the state vector X n . We note that non-standard particle methods designed to minimize the variance of the estimate of p θ (y 0:n ) have recently been proposed [98].
Such quantities are critical when implementing ML parameter estimation procedures; see Section 5. If we substitute p θ ( dx 0:n | y 0:n ) to p θ ( x 0:n | y 0:n ) dx 0:n to approximate S θ n , then we obtain an estimate S θ n which can be computed recursively in time, see e.g. [10,Section 8.3]. For the remainder of this paper we will refer to this approximation as the path space approximation. Even when (3.11) holds, there exists 0 < F θ , G θ < ∞ such that the asymptotic bias [25] and variance [85] satisfy where the variance is w.r.t the law of the particle filter. The fact that the variance grows at least quadratically in time follows from the degeneracy problem and makes S θ n unsuitable for some on-line likelihood based parameter estimation schemes discussed in Section 5.

SMOOTHING
In this section, the parameter θ is still assumed known and we focus on smoothing, that is the problem of estimating the latent variables X 0:T given a fixed batch of observations y 0:T . Smoothing for a fixed parameter θ is at the core of the two main particle ML parameter inference techniques described in Section 5 as these procedures require computing smoothed additive functionals of the form (3.14). Clearly one could unfold the recursion (3.1) from n = 0 to n = T to obtain p θ ( x 0:T | y 0:T ). However, as pointed out in the previous section, the path space approximation (3.7) suffers from the degeneracy problem and yields potentially high variance estimates of (3.14) as (3.15) holds. This has motivated the development of alternative particle approaches to approximate p θ ( x 0:T | y 0:T ) and its marginals.

Fixed-lag approximation
For state-space models with "good" forgetting properties (e.g. (3.11)), we have for L large enough; that is observations collected at times k > n + L do not bring any significant additional information about X 0:n . In particular, when having to evaluate S θ T of the form (3.14) we can approximate the expectation of s n (x n−1:n ) w.r.t p θ ( x n−1:n | y 0:T ) by its expectation w.r.t p θ x n−1:n | y 0:(n+L)∧T .
Algorithmically, a particle implementation of (4.1) means not resampling the components X i 0:n of the particles X i 0:k obtained by particle filtering at times k > n + L. This was first suggested in [59] and used in [10, section 8.3], [78]. This algorithm is simple to implement but the main practical problem is the choice of L. If taken too small, then p θ x 0:n | y 0:(n+L)∧T ) is a poor approximation of p θ ( x 0:n | y 0:T ). If taken too large, the degeneracy remains substantial. Moreover, even as N → ∞, this particle approximation will have a non-vanishing bias since p θ ( x 0:n | y 0:T ) = p θ x 0:n | y 0:(n+L)∧T ) .

Principle
The joint smoothing density p θ ( x 0:T | y 0:T ) can be expressed as a function of the filtering densities {p θ ( x n | y 0:n )} T n=0 using the following key decomposition where p θ ( x n | y 0:n , x n+1 ) is a backward (in time) Markov transition density given by A backward in time recursion for {p θ ( x n | y 0:T )} T n=0 follows by integrating out x 0:n−1 and x n+1:T in (4.2) while applying (4.3) This is referred to as forward-backward smoothing as a forward pass yields {p θ ( x n | y 0:n )} T n=0 which can be used in a backward pass to obtain {p θ ( x n | y 0:T )} T n=0 . Combined to {p θ ( x n | y 0:n , x n+1 )} T −1 n=0 , this allows us to obtain S θ T . An alternative to these forward-backward procedures is the generalized two-filter formula [8].
This Forward Filtering Backward Sampling (FFBSa) procedure was proposed in [48]. It requires O (N (T + 1)) operations to generate a single path X 0:T as sampling from (4.5) costs O (N ) operations. However, as noted in [29], it is possible to sample using rejection from an alternative approximation of p θ ( x n | y 0:n , X n+1 ) in O (1) operations if we use an unweighted particle approximation of p θ ( x n | y 0:n ) in (4.3) and if the transition prior satisfies f θ (x ′ |x) ≤ C < ∞. Hence, with this approach, sampling a path X 0:T costs only on average O (T + 1) operations. A related rejection technique was proposed in [51]. In practice, one may generate N such trajectories to compute Monte Carlo averages that approximates smoothing expectations E [ϕ (X 0:T ) |y 0:T ]. In that scenario, the first approach costs O N 2 (T + 1) , while the second approach costs O (N (T + 1)) on average. In some applications, the rejection sampling procedure can be computationally costly as the acceptance probability can be very small for some particles; see for example Section 4.3 in [79] for empirical results. This has motivated the development of hybrid procedures combining FFBSa and rejection sampling [90]. We can also directly approximate the marginals {p θ ( x n | y 0:T )} T n=0 . Assuming we have an ap- T then by using (4.4) and (4.5), we obtain the approximation p θ ( dx n | y 0: This Forward Filtering Backward Smoothing (FFBSm, where 'm' stands for 'marginal') procedure requires O N 2 (T + 1) operations to approximate {p θ ( x n | y 0:T )} T n=0 instead of O (N (T + 1)) for the path space and fixed-lag methods. However this high computational complexity of forwardbackward estimates can be reduced using fast computational methods [61]. Particle approximations of generalized two-filter smoothing procedures have also been proposed in [8], [39].

Principle
Whenever we are interested in computing the sequence S θ n n≥0 recursively in time, the forward-backward procedure described above is cumbersome as it requires performing a new backward pass with n + 1 steps at time n. An important but not well-known result is that it is possible to implement exactly the forward-backward procedure using only a forward procedure. This result is at the core of [35] but its exposition relies on tools which are non-standard for statisticians. We follow here the simpler derivation proposed in [26] which simply consists of rewriting (3.14) as It can be easily checked using (4.2) that V θ n (x n ) satisfies the following forward recursion for n ≥ 0 . In practice, we shall approximate the function V θ n on a certain grid of values x n , as explained in the next section.

Particle implementation
We can easily provide a particle approximation of the forward smoothing recursion. Assume you have access to approximations V θ n X i . Then when updating our particle filter to obtain , we can directly compute the particle approximations V θ n+1 X i n+1 by plugging (4.5) and p θ ( dx n | y 0:n ) in (4.7)-(4.9) to obtain This approach requires O N 2 (n + 1) operations to compute S θ n at iteration n. A variation over this idea recently proposed in [79] and [93] consists of approximating V θ n+1 X i n+1 by sampling X i,j n ∼ p θ dx n | y 0:n , X i n+1 for j = 1, ..., K to obtain When it is possible to sample from p θ dx n | y 0:n , X i n+1 in O (1) operations using rejection sampling, (4.12) provides a Monte Carlo approximation to (4.10) of overall complexity O (N K).

Convergence results for particle smoothing
Empirically, for a fixed number of particles, these smoothing procedures perform significantly much better than the naive path space approach to smoothing (that is, simply propagating forward the complete state trajectory within a particle filtering algorithm). Many theoretical results validating these empirical findings have been established under assumption (3.11) and additional regularity assumptions. The particle estimate of S θ n based on the fixed-lag approximation (4.1) has an asymptotic variance in n/N with a non-vanishing (as N → ∞) bias proportional to n and a constant decreasing exponentially fast with L [78]. In [26], [27], [29], it is shown that when (3.11) holds, there exists 0 < F θ , H θ < ∞ such that the asymptotic bias and variance of the particle estimate of S θ n computed using the forward-backward procedures satisfy The bias for the path space and forward-backward estimators of S θ n are actually equal [27]. Recently, it has also been established in [79] that, under similar regularity assumptions, the estimate (4.12) also admits an asymptotic variance in n/N whenever K ≥ 2.

MAXIMUM LIKELIHOOD PARAMETER ESTIMATION
We describe in this Section how the particle filtering and smoothing techniques introduced in Sections 3 and 4 can be used to implement maximum likelihood parameter estimation techniques.

Off-line methods
We recall that ℓ T (θ) denote the log-likelihood function associated to data y 0:T introduced in Section 2. So as to maximize ℓ T (θ) , one can rely on standard nonlinear optimization methods, e.g. using quasi-Newton or gradient-ascent techniques. We will limit ourselves to these approaches even if they are sensitive to initialization and might get trapped in a local maximum.

Likelihood function evaluation
We have seen in Section 3 that ℓ T (θ) can be approximated using particle methods, for any fixed θ ∈ Θ. One may wish then to treat ML estimation as an optimization problem using Monte Carlo evaluations of ℓ T (θ). When optimizing a function calculated with a Monte Carlo error, a popular strategy is to make the evaluated function continuous by using common random numbers over different evaluations to ease the optimization. Unfortunately, this strategy is not helpful in the particle context. Indeed, in the resampling stage, particles which admits a piecewise constant and hence discontinuous cumulative distribution function (cdf). A small change in θ will cause a small change in the importance weights {W i n+1 } N i=1 and this will potentially generate a different set of resampled particles. As a result, the log-likelihood function estimate will not be continuous in θ even if ℓ T (θ) is continuous.
To bypass this problem, an importance sampling method was introduced in [52] but it has computational complexity O N 2 (T + 1) and only provides low variance estimates in the neighborhood of a suitably preselected parameter value. In the restricted scenario where X ⊆ R, an elegant solution to the discontinuity problem was proposed in [76]. The method uses common random numbers and introduces a "continuous" version of the resampling step by finding a permutation σ such that X and defining a piecewise linear approximation of the resulting cdf from which particles are resampled, i.e.
This method requires O (N (T + 1) log N ) operations due to the sorting of the particles but the resulting continuous estimate of ℓ T (θ) can be maximized using standard optimization techniques. Extensions to the multivariate case where X ⊆ R nx (with n x > 1) have been proposed in [63] and [23]. However, the scheme [63] does not guarantee continuity of the likelihood function estimate and only provides log-likelihood estimates which are positively correlated for neighboring values in the parameter space whereas the scheme in [23] has O N 2 computational complexity and relies on a non-standard particle filtering scheme.
When θ is high-dimensional, the optimization over the parameter space may be made more efficient if provided with estimates of the gradient. This is exploited by the algorithms described in the forthcoming sections.

Gradient ascent
The log-likelihood ℓ T (θ) may be maximized with the following steepest ascent algorithm: at iteration k + 1 where ∇ θ ℓ T (θ)| θ=θ k is the gradient of ℓ T (θ) w.r.t θ evaluated at θ = θ k and {γ k } is a sequence of positive real numbers, called the step-size sequence. Typically, γ k is determined adaptively at iteration k using a line search or the popular Barzilai-Borwein alternative. Both schemes guarantee convergence to a local maximum under weak regularity assumptions; see [100] for a survey.
The score vector ∇ θ ℓ T (θ) can be computed by using Fisher's identity given in (2.4). Given (2.2), it is easy to check that the score is of the form (3.14). An alternative to Fisher's identity to compute the score is presented in [21] but this also requires computing an expectation of the form (3.14).
These score estimation methods are not applicable in complex scenarios where it is possible to sample from f θ (x ′ |x) but the analytical expression of this transition kernel is unavailable [53]. For those models, a naive approach is to use a finite difference estimate of the gradient; however this might generate too high a variance estimate. An interesting alternative presented in [54], under the name of iterated filtering, consists of deriving an approximation of ∇ θ ℓ T (θ)| θ=θ k based on the posterior moments {E(ϑ n |y 0:n ), V(ϑ n |y 0:n )} T n=0 of an artificial state-space model with latent Markov process {Z n = (X n , ϑ n )} T n=0 , It is shown in [54] that this approximation improves as σ 2 , τ 2 → 0 and σ 2 /τ 2 → 0. Clearly as the variance σ 2 of the artificial dynamic noise {ε n } on the θ−component decreases, it will be necessary to use more particles to approximate ∇ θ l T (θ)| θ=θ k as the mixing properties of the artificial dynamic model deteriorates.

Expectation-Maximization
Gradient ascent algorithms can be numerically unstable as they require to scale carefully the components of the score vector. The Expectation Maximization (EM) algorithm is a very popular alternative procedure for maximizing ℓ T (θ) [24]. At iteration k + 1, we set The sequence {ℓ T (θ k )} k≥0 generated by this algorithm is non-decreasing. The EM is usually favored by practitioners whenever it is applicable as it is numerically more stable than gradient techniques. In terms of implementation, the EM consists of computing a n s -dimensional summary statistic of the form (3.14) when p θ (x 0:T , y 0:T ) belongs to the exponential family, and the maximizing argument of Q(θ k , θ) can be characterized explicitly through a suitable function Λ : R ns → Θ, i.e.

Discussion of particle implementations
The path space approximation (3.7) can be used to approximate the score (2.4) and the summary statistics of the EM algorithm at the computational cost of O(N (T + 1)); see [10, section 8.3], [78], [85]. Experimentally the variance of the associated estimates increases typically quadratically with T [85]. To obtain estimates whose variance increases only typically linearly with T with similar computational cost, one can use the fixed-lag approximation presented in Section 4.1 or a more recent alternative where the path space method is used but the additive functional of interest, which is a sum of terms over n = 0, . . . , T , is approximated by a sum of similar terms which are now exponentially weighted w.r.t n [77]. These methods introduce a non-vanishing asymptotic bias difficult to quantify but appear to perform well in practice.
To improve over the path space method without introducing any such asymptotic bias, the FFBSm and forward smoothing discussed in Sections 4.2 and 4.3 as well as the generalized twofilter smoother have been used [87], [26], [27], [85], [8]. Experimentally the variance of the associated estimates increases typically linearly with T [85] in agreement with the theoretical results in [26], [27], [29]. However the computational complexity of these techniques is O (N 2 (T + 1)). For a fixed computational complexity of order O N 2 (T + 1) , an informal comparison of the performance of the path space estimate using N 2 particles and the forward-backward estimate using N particles suggest that both estimates admit a Mean Square Error (MSE) of order O(N −2 (T + 1)), but the MSE of the path space estimate is variance dominated whereas the forward-backward estimates are bias dominated. This can be understood by decomposing the MSE as the sum of the squared bias and the variance and then substituting appropriately for N 2 particles in (3.15) for the path space method and for N particles in (4.13) for the forward-backward estimates. We confirm experimentally this fact in Section 7.1.
These experimental results suggest that these particle smoothing estimates might thus of limited interest compared to the path based estimates for ML parameter inference when accounting for computational complexity. However, this comparison ignores that the O(N 2 ) computational complexity of these particle smoothing estimates can be reduced to O (N ) by sampling approximately from p θ (x 0:T |y 0:T ) with the FFBSa procedure in Section 4.2 or by using fast computational methods [61]. Related O (N ) approaches have been developed for generalized two-filter smoothing [7], [39]. When applicable, these fast computational methods should be favored.

On-line methods
For a long observation sequence the computation of the gradient of ℓ T (θ) can be prohibitive, moreover we might have real-time constraints. An alternative would be a recursive procedure in which the data is run through once sequentially. If θ n is the estimate of the model parameter after the first n observations, a recursive method would update the estimate to θ n+1 after receiving the new data y n . Several on-line variants of the ML procedures described earlier are now presented. For these methods to be justified, it is crucial for the observation process to be ergodic for the limiting averaged likelihood function ℓ T (θ) /T to have a well-defined limit ℓ (θ) as T → +∞.

On-line gradient ascent
An alternative to gradient ascent is the following parameter update scheme at time n ≥ 0 (5.6) θ n+1 = θ n + γ n+1 ∇ log p θ (y n |y 0:n−1 )| θ=θn , where the positive non-increasing step-size sequence {γ n } n≥1 satisfies n γ n = ∞ and n γ 2 n < ∞ [6], [68]; e.g. γ n = n −α for 0.5 < α ≤ 1. Upon receiving y n , the parameter estimate is updated in the direction of ascent of the conditional density of this new observation. In other words, one recognizes in (5.6) the update of the gradient ascent algorithm (5.1), except that the partial (up 13 to time n) likelihood is used. The algorithm in the present form is however not suitable for online implementation, because evaluating the gradient of log p θ (y n |y 0:n−1 ) at the current parameter estimate requires computing the filter from time 0 to time n using the current parameter value θ n . An algorithm bypassing this problem has been proposed in the literature for a finite state-space latent process in [68]. It relies on the following update scheme (5.7) θ n+1 = θ n + γ n+1 ∇ log p θ 0:n (y n |y 0:n−1 ) where ∇ log p θ 0:n (y n |y 0:n−1 ) is defined as (5.8) ∇ log p θ 0:n (y n |y 0:n−1 ) = ∇ log p θ 0:n (y 0:n ) − ∇ log p θ 0:n−1 (y 0:n−1 ), with the notation ∇ log p θ 0:n (y 0:n ) corresponding to a 'time-varying' score which is computed with a filter using the parameter θ p at time p. The update rule (5.7) can be thought of as an approximation to the update rule (5.6). If we use Fisher's identity to compute this 'time-varying' score, then we have for 1 ≤ p ≤ n The asymptotic properties of the recursion (5.7) (i.e. the behavior of θ n in the limit as n goes to infinity) has been studied in [68] for a finite state-space HMM. It is shown that under regularity conditions this algorithm converges towards a local maximum of the average log-likelihood ℓ (θ), ℓ (θ) being maximized at the 'true' parameter value under identifiability assumptions. Similar results hold for the recursion (5.6).

On-line Expectation-Maximization
It is also possible to propose an on-line version of the EM algorithm. This was originally proposed for finite state-space and linear Gaussian models in [36], [44]; see [13] for a detailed presentation in the finite state-space case. Assume that p θ (x 0:n , y 0:n ) is in the exponential family. In the on-line implementation of EM, running averages of the sufficient statistics n −1 S θ n are computed [12], [36]. Let {θ p } 0≤p≤n be the sequence of parameter estimates of the on-line EM algorithm computed sequentially based on y 0:n−1 . When y n is received, we compute (5.10) S θ 0:n = γ n+1´sn (x n−1:n ) p θ 0:n (x n−1 , x n |y 0:n )dx n−1:n where {γ n } n≥1 needs to satisfy n γ n = ∞ and n γ 2 n < ∞. Then the standard maximization step (5.5) is used as in the batch version The recursive calculation of S θ 0:n is achieved by setting V θ 0 = 0, then computing (5.12) and finally (5.13) S θ 0:n =ˆV θ 0:n (x n ) p θ 0:n (x n |y 0:n )dx n .
Again, the subscript θ 0:n on p θ 0:n (x 0:n |y 0:n ) indicates that the posterior density is being computed sequentially using the parameter θ p at time p ≤ n. The filtering density then is advanced from time n − 1 to time n by using f θn (x n |x n−1 ), g θn (y n |x n ) and p θn (y n |y 0:n ) in the fraction of the rhs of (3.1). Whereas the convergence of the EM algorithm towards a local maximum of the average log-likelihood ℓ (θ) has been established for i.i.d. data [11], its convergence for state-space models remains an open problem despites empirical evidence it does [12], [13], [27]. This has motivated the development of modified versions of the on-line EM algorithm for which convergence results are easier to establish [4], [66]. However, the on-line EM presented here performs empirically usually better [67].

Discussion of particle implementations
Both the on-line gradient and EM procedures require approximating terms (5.8) and (5.10) of the form (3.14), except that the expectation is now w.r.t the posterior density p θ 0:n (x 0:n |y 0:n ) which is updated using the parameter θ p at time p ≤ n. In this on-line framework, only the path space, fixed lag smoothing and forward smoothing estimates are applicable, the fixed lag approximation is also applicable but introduces a non-vanishing bias. For the on-line EM algorithm, similarly to the batch case discussed in Section 5.1.4, the benefits of using the forward smoothing estimate [27] compared to the path space estimate [12] with N 2 particles are rather limited, as experimentally demonstrated in Section 7.1. However for the on-line gradient ascent algorithm, the gradient term ∇ log p θ 0:n (y n |y 0:n−1 ) in (5.7) is a difference between two score-like vectors (5.8) and the behavior of its particle estimates differs significantly from its EM counterpart. Indeed the variance of the particle path estimate of ∇ log p θ 0:n (y n |y 0:n−1 ) increases linearly with n, yielding an unreliable gradient ascent procedure, whereas the particle forward smoothing estimate has a variance uniformly bounded in time under appropriate regularity assumptions and yields a stable gradient ascent procedure [28]. Hence the use of a procedure of computational complexity O N 2 is clearly justified in this context. The very recent paper [93] reports that the computationally cheaper estimate (4.12) appears to exhibit similar properties whenever K ≥ 2 and might prove an attractive alternative.

BAYESIAN PARAMETER ESTIMATION
In the Bayesian setting, we assign a suitable prior density p (θ) for θ and inference is based on the joint posterior density p ( x 0:T , θ| y 0:T ) in the off-line case, or the sequence of posterior densities {p ( x 0:n , θ| y 0:n )} n≥0 in the on-line case.

Particle Markov chain Monte Carlo methods
Using MCMC is a standard approach to approximate p ( x 0:T , θ| y 0:T ). Unfortunately designing efficient MCMC sampling algorithms for nonlinear non-Gaussian state-space models is a difficult task: one-variable-at-a-time Gibbs sampling typically mixes very poorly for such models, whereas blocking strategies that have been proposed in the literature are typically very model-dependent; see for instance [56].
Particle MCMC are a class of MCMC techniques which rely on particle methods to build efficient high dimensional proposal distributions in a generic manner [3]. We limit ourselves here to the presentation of the Particle Marginal Metropolis-Hastings (PMMH) sampler, which is an approximation of an ideal MMH sampler for sampling from p ( x 0:T , θ| y 0:T ) which would utilize the following proposal density where q ( θ ′ | θ) is a proposal density to obtain a candidate θ ′ when we are at location θ. The acceptance probability of this sampler is Unfortunately this ideal algorithm cannot be implemented as we cannot sample exactly from p θ ′ ( x 0:T | y 0:T ) and we cannot compute the likelihood terms p θ (y 0:T ) and p θ ′ (y 0:T ) appearing in the acceptance probability. The PMMH sampler is an approximation of this ideal MMH sampler which relies on the particle approximations of these unknown terms. Given θ and a particle approximation p θ (y 0:T ) of p θ (y 0:T ), we sample θ ′ ∼ q ( θ ′ | θ) then run a particle filter to obtain approximations p θ ′ ( dx 0:T | y 0:T ) and p θ ′ (y 0:T ) of p θ ′ ( dx 0:T | y 0:T ) and p θ ′ (y 0:T ). We then sample X ′ 0:T ∼ p θ ′ ( dx 0:T | y 0:T ), that is we choose randomly one of N particles generated by the particle filter, with probability W i T for particle i, and accept (θ ′ , X ′ 0:T ) (and p θ ′ (y 0:T )) with probability The acceptance probability (6.3) is a simple approximation of the "ideal" acceptance probability (6.2). This algorithm was first proposed as an heuristic to sample from p ( θ| y 0:T ) in [41]. Its remarkable feature established in [3] is that it does admit p ( x 0:T , θ| y 0:T ) as invariant distribution whatever being the number of particles N used in the particle approximation [3]. However the choice of N has an impact on the performance of the algorithm. Using large values of N usually results in PMMH averages with variances lower than the corresponding averages using fewer samples but the computational cost of constructing p θ (y 0:T ) increases with N . A simplified analysis of this algorithm suggests that N should be selected such that the standard deviation of the logarithm of the particle likelihood estimate should be around 0.9 if the ideal MMH sampler was using the perfect proposal q(θ ′ |θ) = p (θ ′ |y 0:n ) [83] and around 1.8 if one uses an isotropic normal random walk proposal, the target is a product of d i.i.d. components and d → ∞ [88]. For general proposal and target densities, a recent theoretical analysis and empirical results suggest that this standard deviation should be selected around 1.2 − 1.3 [34]. As the variance of this estimate typically increases linearly with T , this means that the computational complexity is of order O(T 2 ) by iteration.
A particle version of the Gibbs sampler is also available [3] which mimicks the two-component Gibbs sampler sampling iteratively from p ( θ| x 0:T , y 0:T ) and p θ ( x 0:T | y 0:T ). These algorithms rely on a non-standard version of the particle filter where N −1 particles are generated conditional upon a "fixed" particle. Recent improvements over this particle Gibbs sampler introduce mechanisms to rejuvenate the fixed particle, using forward or backward sampling procedures [95], [70], [96]. These methods perform empirically extremely well but, contrary to the PMMH, it is still unclear how one should scale N with T .

On-line Methods
In this context, we are interested in approximating on-line the sequence of posterior densities {p ( x 0:n , θ| y 0:n )} n≥0 . We emphasize that, contrary to the on-line ML parameter estimation procedures, none of the methods presented in this section bypass the particle degeneracy problem. This should come as no surprise. As discussed in Section 3.2.2, even for a fixed θ, the particle estimate of p θ (y 0:n ) has a relative variance that increases linearly with n under favorable mixing assumptions. The methods in this section attempt to approximate p (θ|y 0:n ) ∝ p θ (y 0:n )p(θ). This is a harder problem as it implicitly requires having to approximate p θ i (y 0:n ) for all the particles θ i approximating p (θ|y 0:n ).

Augmenting the state with the parameter
At first sight, it seems that estimating the sequence of posterior densities {p ( x 0:n , θ| y 0:n )} n≥0 can be easily achieved using standard particle methods, by merely introducing the extended state Z n = (X n , θ n ), with initial density p (θ 0 ) µ θ 0 (x 0 ) and transition density f θn ( x n | x n−1 ) δ θ n−1 (θ n ) ; i.e. θ n = θ n−1 . However, this extended process Z n clearly does not possess any forgetting property (as discussed in Section 3), so the algorithm is bound to degenerate. Specifically, the parameter space is explored only in the initial step of the algorithm. Then, each successive resampling step reduces the diversity of the sample of θ values; after a certain time n, the approximation p ( dθ| y 0:n ) contains a single unique value for θ. This is clearly a poor approach. Even in the much simpler case when there is no latent variable X 0:n , it is shown in [18,Theorem 4] that the asymptotic variance of the corresponding particle estimates diverges at least at a polynomial rate, which grows with the dimension of θ.
A pragmatic approach that has proven useful in some applications is to introduce artificial dynamics for the parameter θ [58] where {ε n } n≥0 is an artificial dynamic noise with decreasing variance. Standard particle methods can now be applied to approximate {p ( x 0:n , θ 0:n | y 0:n )} n≥0 . A related kernel density estimation method also appeared in [71], which proposes to use a kernel density estimate p (θ|y 0:n ) from which one samples from. As before the static parameter is transformed to a slowly time-varying one, whose dynamics is related to the kernel bandwidth. To mitigate the artificial variance inflation, a shrinkage correction is introduced. An improved version of this method has been recently proposed in [42]. It is difficult to quantify how much bias is introduced in the resulting estimates by the introduction of this artificial dynamics. Additionally, these methods require a significant amount of tuning, e.g. choosing the variance of the artificial dynamic noise or the kernel width. However they can perform satisfactorily in practice [42], [71].

Practical filtering
The practical filtering approach proposed in [84] relies on the following fixed-lag approximation (6.5) p ( x 0:n−L , θ| y 0:n−1 ) ≈ p ( x 0:n−L , θ| y 0:n ) for L large enough; that is observations coming after n − 1 presumably brings little information on x 0:n−L . To sample approximately from p ( θ| y 0:n ), one uses the following iterative process: at time n, several MCMC chains are run in parallel to sample from p x n−L+1:n , θ| y 0:n , X i 0:n−L = p x n−L+1:n , θ| y n−L+1:n , X i n−L where the X i n−L have been obtained at the previous iteration, and are such that (approximately) X i n−L ∼ p(x n−L |y 0:n−1 ) ≈ p(x n−L |y 0:n ). Then one collects the first component X i n−L+1 of the simulated sample X i n−L+1:n , increments the time index and runs several new MCMC chains in parallel to sample from p x n−L+2:n+1 , θ| y n−L+2:n+1 , X i n−L+1 and so on. The algorithm is started at time L − 1, with MCMC chains that target p(x 0:L−1 |y 0:L−1 ). Like all methods based on fixed-lag approximation, the choice of the lag L is difficult and this introduces a non-vanishing bias which is difficult to quantify. However, the method performs well on the examples presented in [84].

Using MCMC steps within particle methods
To avoid the introduction of an artificial dynamic model or of a fixed-lag approximation, an approach originally proposed independently in [37] and [47] consists of adding MCMC steps to re-introduce "diversity" among the particles. Assume we use an auxiliary particle filter to approximate {p ( x 0:n , θ| y 0:n )} n≥0 then the particles X i 0:n , θ i n obtained after the sampling step at time n are approximately distributed according to p ( x 0:n , θ| y 0:n ) ∝ p ( x 0:n−1 , θ| y 0:n−1 ) q θ (x n , y n |x n−1 ) .
If we use the SISR algorithm, then we can alternatively use an MCMC step of invariant density p ( x 0:n , θ| y 0:n ) after the resampling step at time n. Contrary to standard applications of MCMC, the kernel does not have to be ergodic. Ensuring ergodicity would indeed require one to sample an increasing number of variables as n increasesthis algorithm would have an increasing cost per iteration, which would prevents its use in on-line scenarios but it can be an interesting alternative to standard MCMC and was suggested in [65]. In practice one therefore sets X i 0:n−L = X i 0:n−L and only sample θ i and X i n−L+1:n , where L is a small integer; often L = 0 (only θ is updated). Note that the memory requirements for this method do not increase over time if p θ (x 0:n , y 0:n ) is in the exponential family and thus can be summarized by a set of fixed dimensional sufficient statistics s n (x 0:n , y 0:n ). This type of methods was first used in to perform on-line Bayesian parameter estimation in a context where p θ (x 0:n , y 0:n ) is in the exponential family [47], [37]. Similar strategies were adopted in [2] and [89]. In the particular scenario where q θ (x n |y n , x n−1 ) = p θ (x n |y n , x n−1 ) and q θ (y n |x n−1 ) = p θ (y n |x n−1 ), this method was mentioned in [2], [91] and is discussed at length in [74] who named it particle learning. Extensions of this strategy to parameter estimation in conditionally linear Gaussian models, where a part of the state is integrated out using Kalman techniques [16], [31], is proposed in [14].
As opposed to the methods relying on kernel or artificial dynamics, these MCMC-based approaches have the advantage of adding diversity to the particles approximating p (θ|y 0:n ) without perturbing the target distribution. Unfortunately, these algorithms rely implicitly on the particle approximation of the density p (x 0:n |y 0:n ) even if algorithmically it is only necessary to store some fixed-dimensional sufficient statistics s n (X i 0:n , y 0:n ) . Hence in this respect they suffer from the degeneracy problem. This was noticed as early as in [2]; see also the word of caution in the conclusion of [37], [4] and [20]. The practical implications are that one observes empirically that the resulting Monte Carlo estimates can display quite a lot of variability over multiple runs as demonstrated in Section 7.2. This should not come as a surprise as the sequence of posterior distributions does not have exponential forgetting properties, hence there is an accumulation of Monte Carlo errors over time.

The SMC 2 algorithm
The SMC 2 algorithm introduced simultaneously in [19] and [45] may be considered as the particle equivalent of Particle MCMC. It mimics an "ideal" particle algorithm proposed in [17] approximating sequentially {p ( θ| y 0:n )} n≥0 where N θ particles (in the θ-space) are used to explore these distributions. The N θ particles at time n are reweighted according to p θ (y 0:n+1 )/p θ (y 0:n ) at time n + 1. As these likelihood terms are unknown, we substitute to them p θ (y 0:n+1 )/p θ (y 0:n ) wherep θ (y 0:n ) is a particle approximation of the partial likelihood p θ (y 0:n ), obtained by a running a particle filter of N x particles in the x−dimension, up to time n, for each 18 of the N θ θ−particles. When particle degeneracy (in the θ−dimension) reaches a certain threshold, θ−particles are refreshed through the succession of a resampling step, and an MCMC step, which in these particular settings takes the form of a PMCMC update. The cost per iteration of this algorithm is not constant and, additionally, it is advised to increase N x with n for the relative variance ofp θ (y 0:n ) not to increase, therefore it cannot be used in truly on-line scenarios. Yet there are practical situations where it may be useful to approximate jointly all the posteriors p(θ|y 1:n ), for 1 ≤ n ≤ T , for instance to assess the predictive power of the model.

EXPERIMENTAL RESULTS
We focus on illustrating numerically a few algorithms and the impact of the degeneracy problem on parameter inference. This last point is motivated by the fact that particle degeneracy seems to have been overlooked by many practitioners. In this way numerical results may provide valuable insights.
We will consider the following simple scalar linear Gaussian state space model: where V n , W n are independent zero-mean and unit-variance Gaussians and ρ ∈ [−1, 1]. The main reason for choosing this model is that Kalman filter recursions can be implemented to provide the exact values of the summary statistics S θ n used for ML estimation through the EM algorithm, and to compute the exact likelihood p θ (y 0:n ). Hence, using a fine discretization of the low-dimensional parameter space, we can compute a very good approximation of the true posterior density p ( θ| y 0:n ). In this model it is straightforward to present numerical evidence of some effects of degeneracy for parameter estimation and show how it can be overcome by choosing an appropriate particle method.

Maximum likelihood methods
As ML methods require approximating smoothed additive functionals S θ n of the form (3.14), we begin by investigating the empirical bias, variance and MSE of two standard particle estimates of S θ n , where we set s k (x k−1 , x k ) = x k−1 x k for the model described in (7.1). The first estimate relies on the path space method with computational cost O(N ) per time, which uses p θ ( dx 0:n | y 0:n ) in (3.7) to approximate S θ n as S θ n ; see [10,Section 8.3] for more details. The second estimate relies on the forward implementation of FFBSm presented in Section 4.3 using (4.7)-(4.11); see [27]. Recall that this procedure has a computational cost that is O(N 2 ) per time for N particles and provides the same estimates as the standard forward-backward implementation of FFBSm. For the sake of brevity we will not consider the remaining smoothing methods of Section 4; for the fixed-lag and the exponentially weighted approximations we refer the reader to [78] respectively [77] for numerical experiments.
We use a simulated dataset of size 6 × 10 4 obtained using θ * = (ρ * , τ 2 * , σ 2 * ) = (0.8, 0.1, 1) and then generate 300 independent replications of each method in order to compute the empirical bias and variance ofŜ θ * n when θ is fixed to θ * . In order to make a comparison that takes into account the computational cost, we use N 2 particles for the O(N ) method and N for the O(N 2 ) one. We look separately at the behavior of the bias ofŜ θ n and the variance and MSE of the rescaled estimateŝ S θ n / √ n. The results are presented in Figure 1 for N = 50, 100, 200.
For both methods the bias grows linearly with time, this growth being higher for the O(N 2 ) method. For the variance ofŜ θ n / √ n, we observe a linear growth with time for the O(N ) method with N 2 particles whereas this variance appears roughly constant for the O(N 2 ) method. Finally, the MSE ofŜ θ n / √ n grows for both methods linearly as expected. In this particular scenario, the constants of proportionality are such that the MSE is lower for the O(N ) method than for the O(N 2 ) method. In general, we can expect that the O(N ) method be superior in terms of the bias and the O(N 2 ) method superior in terms of the variance. These results are in agreement with the theoretical results in the literature [26], [27], [29] but additionally show that the lower bound on the variance growth ofŜ θ n for the O(N ) method of [85] appears sharp. We proceed to see how the bias and variance of the estimates of S θ n affect the ML estimates, when the former are used within both an off-line and an on-line EM algorithm; see Figures 2 and 3 respectively. For the model in (7.1) the E-step corresponds to computing S θ n where s k (x k−1 , x k ) = (y k − x k ) 2 , x 2 k−1 , x k−1 x k , x 2 k and the M-step update function is given by We compare the estimates of θ * when the E-step is computed using the O(N ) and the O(N 2 ) methods described in the previous section with 150 2 and 150 particles respectively. A simulated dataset for θ * = (ρ * , τ * , σ * ) = (0.8, 1, .2) will be used. In every case we will initialize the algorithm using θ 0 = (0.1, 0.1, 0.2) and assume σ * is known. In Figures 2 and 3 we present the results obtained using 150 independent replications of the algorithm. For the off-line EM, we use 25 iterations for T = 100, 1000, 2500, 5000, 10000. For the on-line EM, we use T = 10 5 with the step size set as γ n = n −0.8 and for the first 50 iterations no M-step update is performed. This "freezing" phase is required to allow for a reasonable estimation of the summary statistic; see [12], [13] for more details. Note that in Figure 3 we plot only the results after the algorithm has converged, i.e. for 21 n ≥ 5×10 4 . In each case, both the O(N ) and the O(N 2 ) methods yield fairly accurate results given the low number of particles used. However we note, as observed previously in the literature, that the on-line EM as well as the on-line gradient ascent method requires a substantial number of observations, i.e. over 10000, before achieving convergence [12], [13], [27], [85]. For smaller datasets, these algorithms can also be used by going through the data say K times. Typically this method is cheaper than iterating (5.1) or (5.4)-(5.5) K times the off-line algorithms and can yield comparable parameter estimates [99]. Experimentally, the properties of the estimates of S θ n discussed earlier appear to translate into properties of the resulting parameter estimates: the O(N ) method provides estimates with less bias but more variance than the O(N 2 ) method.
For more numerical examples regarding the remaining methods discussed in Section 5, we refer the reader to [53], [54] for iterated filtering, to [26], [27], [85] for comparisons of the O(N ) and O(N 2 ) methods for EM and gradient ascent, to [12] for the O(N ) on-line EM, to [76] and [63,Chap. 10] for smooth likelihood function methods and to [10, for a detailed exposition of off-line EM methods.

Bayesian methods
We still consider the model in (7.1) but simplify it further by fixing either ρ or τ . This is done in order to keep the computations of the benchmarks that use Kalman computations on a grid relatively inexpensive. For those parameters that are not fixed, we shall use the following independent priors: a uniform on [−1, 1] for ρ, and inverse gamma for τ 2 , σ 2 with the shape and scale parameter pair being (a, b) and (c, d) respectively with a = b = c = d = 1. In all the subsequent examples, we will initialize the algorithms by sampling θ from the prior.  i.e. empirical variance (over independent runs) for the estimator of the mean of p(τ 2 |y0:n) using particle method with MCMC normalized with the true posterior variance computed using Kalman filtering on a grid. Bottom left: average (over independent runs) of the estimated variance of p(τ 2 |y0:n) using particle method with MCMC normalized with the true posterior variance. Bottom right: variance of the {p (y0:n)} n≥0 ; All plots are computed using N = 5000 and over 100 different independent runs.
We proceed to examine the combination of particle method with MCMC methods described in Section 6.2.3. We focus on an efficient implementation of this idea discussed in [74] which can be put in practice for the simple model under consideration. We investigate the effect of the degeneracy problem in this context. The numerical results obtained in this section have been produced in Matlab (code available from the first author), and double-checked using the R program available on the personal web page of the first author of [74].
We first focus of the estimate of the posterior of θ = (τ 2 , σ 2 ) given a long sequence of simulated observations with τ = σ = 1. In this scenario, p θ (x 0:n , y 0:n ) admits the following two-dimensional sufficient statistics, s n (x 0:n , y 0:n ) = n k=1 (x k − x k−1 ) 2 , n k=0 (y k − x k ) 2 , and θ can be updated using Gibbs steps. We use T = 5 × 10 4 and N = 5000. We ran the algorithm over 100 independent runs over the same dataset. We present the results only for τ 2 and omit the ones for σ 2 as these were very similar. The top left panel of Figure 4 shows the box plots for the estimates of the posterior mean, and the top right panel shows how the corresponding relative variance of the estimator for the posterior mean evolves with time. Here the relative variance is defined as the ratio of the empirical variance (over different independent runs) of the posterior mean estimates at time n over the true posterior variance at time n, which in this case is approximated using a Kalman filter on a fine grid. This quantity exhibits a steep increasing trend when n ≥ 15000 and confirms the aforementioned variability of the estimates of the posterior mean. In the bottom left panel of Figure 4 we plot the average (over different runs) of the estimators of the variance of p(θ|y 0:n ). This average variance is also scaled/normalized by the actual posterior variance. The latter is again computed using Kalman filtering on a grid. This ratio between the average estimated variance of the posterior over the true one decreases with time n and it shows that the supports of the approximate posterior densities provided by this method cover on average only a small portion of the support of the true posterior. These experiments confirm that in this example the particle method with MCMC steps fails to adequately explore the space of θ. Although the box plots provide some false sense of security, the relative and scaled average variance clearly indicate that any posterior estimates obtained from a single run of particle method with MCMC steps should be used with caution. Furthermore, in the the bottom right panel of Figure 4 we also investigate experimentally the empirical variance of the marginal likelihood estimates {p (y 0:n )} n≥0 . This variance appears to increase quadratically with n for the particle method with MCMC moves instead of linearly as it does for state-space models with good mixing properties. This suggests that to one should increase the number of particles quadratically with the time index to obtain an estimate of the marginal likelihood whose relative variance remains uniformly bounded with respect to the time index. Although we attribute this quadratic variance growth to the degeneracy problem, the estimatep (y 0:n ) is not the particle approximation of a smoothed additive functional, thus there is not yet any theoretical convergence result explaining rigorously this phenomenon.
One might argue that these particle methods with MCMC moves are meant to be used with larger N and/or shorter data sets T . We shall consider this time a slightly different example where τ = 0.1 is known and we are interested in estimating the posterior of θ = (ρ, σ 2 ) given a sequence of observations obtained using ρ = 0.5 and σ = 1. In that case, the sufficient statistics are s n (x 0:n , y 0:n ) = n k=1 x k−1 x k , n−1 k=0 x 2 k−1 , n k=0 (y k − x k ) 2 , and the parameters can be rejuvenated through a single Gibbs update. In addition, we let T = 5000 and use N = 10 4 particles. In Figures 5 we display the estimated marginal posteriors p ( ρ| y 0:n ) and p σ 2 y 0:n obtained from 50 independent replications of the particle method. On this simple problem, the estimated posteriors seem consistently rather inaccurate for ρ, whereas they perform better for σ 2 but with some non-negligible variability over runs which increases as T increases. Similar observations have been reported in [20] and remain unexplained: for some parameters this methodology appears to provide reasonable results despite the degeneracy problem and for others it provides very unreliable results.
We investigate further the performance of this method in this simple example by considering the same example for T = 1000 but now consider two larger numbers of particles, N = 7.5 × 10 4 and N = 6 × 10 5 , over 50 different runs. Additionally we compare the resulting estimates with estimates provided by the particle Gibbs sampler of [70] using the same computational cost, that is N = 50 particles with 3000 and 24000 iterations respectively. The results are displayed in Figure 6 and 7. As expected, we improve the performance of the particle with MCMC moves when N increases for a fixed time horizon T . For a fixed computational complexity, the particle Gibbs sampler estimates appear to display less variability. For a higher dimensional parameter θ and/or very vague priors, this comparison would be more favorable to the particle Gibbs sampler as illustrated in [3, pp. 336-338].

CONCLUSION
Most particle methods proposed originally in the literature to perform inference about static parameters in general state-space models were computationally inefficient as they suffered from the degeneracy problem. Several approaches have been proposed to deal with this problem by either adding an artificial dynamic on the static parameter [43], [58], [71] or introducing a fixed-lag approximation [59], [78], [84]. These methods can work very well in practice but it remains unfortunately difficult/impossible to quantify the bias introduced in most realistic applications. Various asymptotically bias-free methods with good statistical properties and a reasonable computational cost have recently appeared in the literature.
To perform batch ML estimation, the forward filter backward sampler/smoother and generalized two filter procedures are recommended whenever the O(N 2 T ) computational complexity per iteration of their direct implementations can be lowered to O(N T ) using, for example, the methods described in [7], [29], [39], [61]. Otherwise, besides a lowering of memory requirements, not much can be gained from these techniques compared to simply using a standard particle filter with N 2 particles. In an on-line ML context, the situation is markedly different. Whereas for the on-line EM algorithm, the forward smoothing approach in [85], [27] of complexity O(N 2 ) per time step will be similarly of limited interest compared to a standard particle filter using N 2 particles, it is crucial to use this approach when performing on-line gradient ascent as demonstrated empirically and established theoretically in [28]. In on-line scenarios where one can admit a random computational complexity at each time step, the method presented in [79] is an interesting alternative when it is applicable. Empirically these on-line ML methods converge rather slowly and will be primarily be useful for large datasets.
In a Bayesian framework, batch inference can be conducted using particle MCMC methods [3], [70]. However these methods are computationally expensive as, for example, an efficient implementation of the PMMH has a computational complexity of order O(T 2 ) per iteration [34]. On-line Bayesian inference remains a challenging open problem as all methods currently available, including particle methods with MCMC moves [37], [89], [14] suffer from the degeneracy problem. These methods should not be ruled out but should be used cautiously as they can provide unreliable results even in simple scenarios as demonstrated in our experiments.
Very recent papers in this dynamic research area have proposed to combine individual parameter estimation techniques so as to design more efficient inference algorithms. For example, [22] suggests to use the score estimation techniques developed for ML parameter estimation to design better proposal distributions for the PMMH algorithm whereas [40] demonstrates that particle methods with MCMC moves might be fruitfully used in batch scenarios when plugged into a particle MCMC scheme.