Coupling stochastic EM and approximate Bayesian computation for parameter inference in state-space models

We study the class of state-space models and perform maximum likelihood estimation for the model parameters. We consider a stochastic approximation expectation–maximization (SAEM) algorithm to maximize the likelihood function with the novelty of using approximate Bayesian computation (ABC) within SAEM. The task is to provide each iteration of SAEM with a filtered state of the system, and this is achieved using an ABC sampler for the hidden state, based on sequential Monte Carlo methodology. It is shown that the resulting SAEM-ABC algorithm can be calibrated to return accurate inference, and in some situations it can outperform a version of SAEM incorporating the bootstrap filter. Two simulation studies are presented, first a nonlinear Gaussian state-space model then a state-space model having dynamics expressed by a stochastic differential equation. Comparisons with iterated filtering for maximum likelihood inference, and Gibbs sampling and particle marginal methods for Bayesian inference are presented.


Introduction
State-space models [Cappé et al., 2005] are widely applied in many fields, such as biology, chemistry, ecology, etc.Let us now introduce some notation.Consider an observable, discrete-time stochastic process {Y t } t≥t 0 , Y t ∈ Y ⊆ R dy and a latent and unobserved continuous-time stochastic process {X t } t≥t 0 , X t ∈ X ⊆ R dx .Process X t ∼ p(x t |x s , θ x ) is assumed Markov with transition densities p(•), s < t.Processes {X t } and {Y t } depend on their own (unknown) vector-parameters θ x and θ y , respectively.We consider {Y t } as a measurement-error-corrupted version of {X t } and assume that observations for {Y t } are conditionally independent given {X t }.The state-space model can be summarised as where X 0 ≡ X t 0 .We assume f (•) a known density (or probability mass) function set by the modeller.Regarding the transition density p(x t |x s , •), this is typically unknown except for very simple toy models.Goal of our work is to estimate the parameters (θ x , θ y ) by maximum likelihood, using observations Y 1:n = (Y 1 , ..., Y n ) collected at discrete times {t 1 , ..., t n }.Here Y j ≡ Y t j and we use z 1:n to denote a generic sequence (z 1 , ..., z n ).For ease of notation we refer to the vector θ := (θ x , θ y ) as the object of our inference.
Parameters inference for state-space models has been widely developed, and sequential Monte Carlo (SMC) methods are now considered the state-of-art when dealing with nonlinear/non-Gaussian state space models (see Kantas et al., 2015 for a review).Methodological advancements have especially considered Bayesian approaches.In Bayesian inference the goal is to derive analytically the posterior distribution π(θ|Y 1:n ) or, most frequently, implement an algorithm for sampling draws from the posterior.Sampling procedures are often carried out using Markov chain Monte Carlo (MCMC) or SMC embedded in MCMC procedures, see Andrieu and Roberts [2009] and Andrieu et al. [2010].
In this work we instead aim at maximum likelihood estimation of θ.Several methods for maximum likelihood inference in state-space models have been proposed in the literature, including the well-known EM algorithm [Dempster et al., 1977].The EM algorithm computes the conditional expectation of the complete-likelihood for the pair ({Y t }, {X t }) and then produces a (local) maximizer for the likelihood function based on the actual observations Y 1:n .One of the difficulties is how to compute the conditional expectation of the state {X t } given the observations Y 1:n .This conditional expectation can be computed exactly with the Kalman filter when the state-space is linear and Gaussian [Cappé et al., 2005], otherwise it has to be approximated.In this work we focus on a stochastic approximation of {X t }.Therefore, we resort to a stochastic version of the EM algorithm, namely the Stochastic Approximation EM (SAEM) [Delyon et al., 1999].The problem is to generate, conditionally on the current value of θ during the EM maximization, an appropriate "proposal" for the state {X t }, and we use SMC to obtain such proposal.SMC algorithms [Doucet et al., 2001] have already been coupled to stochastic EM algorithms (see e.g.Huys et al. [2006], Huys and Paninski [2009], Lindsten [2013], Ditlevsen and Samson [2014] and references therein).The simplest and most popular SMC algorithm, the bootstrap filter [Gordon et al., 1993], is easy to implement and very general, explicit knowledge of the density f (y t |X t , •) being the only requirement.Therefore the bootstrap filter is often a go-to option for practitioners.Alternatively, in order to select a path {X t } to feed SAEM with, in this paper we follow an approach based on approximate Bayesian computation (ABC) and specifically we use the ABC-SMC method for state-estimation proposed in Jasra et al. [2012].We do not merely consider the algorithm by Jasra et al. within SAEM, but show in detail and discuss how SAEM-ABC-SMC (shortly SAEM-ABC) can in some cases outperform SAEM coupled with the bootstrap filter.
We illustrate our SAEM-ABC approach for (approximate) maximum likelihood estimation of θ using two case studies, a nonlinear Gaussian state-space model and a more complex state-space model based on stochastic differential equations.We also compare our method with the iterated filtering for maximum likelihood estimation [Ionides et al., 2015], Gibbs sampling and particle marginal methods for Bayesian inference [Andrieu and Roberts, 2009] and will also use a special SMC proposal function for the specific case of stochastic differential equations [Golightly and Wilkinson, 2011].
The paper is structured as follows: in section 2 we introduce the standard SAEM algorithm and basic notions of ABC.In section 3 we propose a new method for maximum likelihood estimation by integrating an ABC-SMC algorithm within SAEM.Section 4 shows simulation results and section 5 summarize conclusions.An appendix includes technical details pertaining the simulation studies.Software code can be found online either at https://github.com/umbertopicchini/SAEM-ABC or as supplementary material in the version of this paper published on Computational Statistics (doi:10.1007/s00180-017-0770-y) 1 .

The complete likelihood and stochastic approximation EM
Recall that Y 1:n = (Y 1 , ..., Y n ) denotes the available data collected at times (t 1 , ..., t n ) and denote with X 1:n = (X 1 , ..., X n ) the corresponding unobserved states.We additionally set X 0:n = (X 0 , X 1:n ) for the vector including an initial (fixed or random) state X 0 , that is X 1 is generated as X 1 ∼ p(x 1 |x 0 ).When the transition densities p(x j |x j−1 ) are available in closed form (j = 1, ..., n), the likelihood function for θ can be written as (here we have assumed a random initial state with density p(X 0 )) p(Y 1:n ; θ) = p Y,X (Y 1:n , X 0:n ; θ) dX 0:n = p Y|X (Y 1:n |X 0:n ; θ)p X (X 0:n ; θ) dX 0:n = p(X 0 ) where p Y,X is the "complete data likelihood", p Y|X is the conditional law of Y given X, p X is the law of X, f (Y j |X j ; •) is the conditional density of Y j as in (1) and p X (X 0:n ; •) the joint density of X 0:n .The last equality in (2) exploits the notion of conditional independence of observations given latent states and the Markovian property of {X t }.In general the likelihood (2) is not explicitly known either because the integral is multidimensional and because expressions for transition densities are typically not available except for trivial toy models.
In addition, when an exact simulator for the solution of the dynamical process associated with the Markov process {X t } is unavailable, hence it is not possible to sample from p(X j |X j−1 ; θ), numerical discretization methods are required.Without loss of generality, say that we have equispaced sampling times such that t j = t j−1 + ∆, with ∆ > 0. Now introduce a discretization for the interval [t 1 , t n ] given by {τ 1 , τ h , ..., τ Rh , ..., τ nRh } where h = ∆/R and R ≥ 1.We take τ 1 = t 1 , τ nRh = t n and therefore τ i ∈ {t 1 , ...., t n } for i = 1, Rh, 2Rh, ..., nRh.We denote with N the number of elements in the discretisation {τ 1 , τ h , ..., τ Rh , ..., τ nRh } and with X 1:N = (X τ 1 , . . ., X τ N ) the corresponding values of {X t } obtained when using a given numerical/approximated method of choice.Then the likelihood function becomes where the product in j is over the X t j and the product in i is over the X τ i .

The standard SAEM algorithm
The EM algorithm introduced by Dempster et al. [1977] is a classical approach to estimate parameters by maximum likelihood for models with non-observed or incomplete data.Let us briefly cover the EM principle.The complete data of the model is (Y 1:n , X 0:N ), where X 0:N ≡ X 0:n if numerical discretization is not required, and for ease of writing we denote this as (Y, X) in the remaining of this section.The EM algorithm maximizes the function is the log-likelihood of the complete data and E is the conditional expectation under the conditional distribution p X|Y (•; θ ).At the k-th iteration of a maximum (user defined) number of evaluations , whereas the M-step updates θ(k−1) by maximizing Q k (θ).For cases in which the E-step has no analytic form, Delyon et al. [1999] introduce a stochastic version of the EM algorithm (SAEM) which evaluates the integral Q k (θ) by a stochastic approximation procedure.The authors prove the convergence of this algorithm under general conditions if L c (Y, X; θ) belongs to the regular exponential family where ., . is the scalar product, Λ and Γ are two functions of θ and S c (Y, X) is the minimal sufficient statistic of the complete model.The E-step is then divided into a simulation step (Sstep) of the missing data X (k) under the conditional distribution p X|Y (•; θ(k−1) ) and a stochastic approximation step (SA-step) of the conditional expectation, using (γ k ) k≥1 a sequence of real numbers in [0, 1], such that at each iteration by the value s k defined recursively as follows The M-step is thus the update of the estimates θ(k−1) The starting s 0 can be set to be a vector of zeros.The procedure above can be carried out iteratively for K iterations.The proof of the convergence of SAEM requires the sequence γ k to be such that A typical choice is to consider a warmup period with γ k = 1 for the first K 1 iterations and then Parameter K 1 has to be chosen by the user.However inference results are typically not very sensitive to this tuning parameter.Typical values are K 1 = 250 or 300 and K = 400, see for example Lavielle [2014].Usually, the simulation step of the hidden trajectory X (k) conditionally on the observations Y cannot be directly performed.Lindsten [2013] and Ditlevsen and Samson [2014] use a sequential Monte Carlo algorithm to perform the simulation step for state-space models.We propose to resort to approximate Bayesian computation (ABC) for this simulation step.
It can be noted that the generation of s k in (3) followed by corresponding parameter estimates θ(k) in ( 4) is akin to two steps of a Gibbs sampling algorithm, except that here θ(k) is produced by a deterministic step (for given s k ).We comment further on this aspect in section 4.1.4.

The SAEM algorithm coupled to an ABC simulation step
Approximate Bayesian Computation (ABC, see Marin et al., 2012 for a review) is a class of probabilistic algorithms allowing sampling from an approximation of a posterior distribution.The most typical use of ABC is when posterior inference for θ is the goal of the analysis and the purpose is to sample draws from the approximate posterior π δ (θ|Y).Here and in the following Y ≡ Y 1:n .The parameter δ > 0 is a threshold influencing the quality of the inference, the smaller the δ the more accurate the inference, and π δ (θ|Y) ≡ π(θ|Y) when δ = 0. However in our study we are not interested in conducting Bayesian inference for θ.We will use ABC to sample from an approximation to the posterior distribution π(X 0:N |Y; θ) ≡ p(X 0:N |Y; θ), that is for a fixed value of θ, we wish to sample from π δ (X 0:N |Y; θ) (recall from section 2.1 that when feasible we can take N ≡ n).For simplicity of notation, in the following we avoid specifying the dependence on the current value of θ, which has to be assumed as a deterministic unknown.There are several ways to generate a "candidate" X * 0:N : for example we might consider "blind" forward simulation, meaning that X * 0:N is simulated from p X (X 0:N ) and therefore unconditionally on data (i.e. the simulator is blind with respect to data).Then X * 0:N is accepted if the corresponding Y * simulated from f (•|X * 1:n ) is "close" to Y, according to the threshold δ, where X * 1:n contains the interpolated values of X * 0:N at sampling times {t 1 , ..., t n } and Y * ≡ Y * 1:n .Notice that the appeal of the methodology is that knowledge of the probabilistic features of the data generating model is not necessary, meaning that even if the transition densities p(X j |X j−1 ) are not known (hence p X is unknown) it is enough to be able to simulate from the model (using a numerical scheme if necessary) hence draws X * 0:N are produced by forward-simulation regardless the explicit knowledge of the underlying densities.
Algorithm 1 illustrates a generic iteration k of a SAEM-ABC method, where the current value of the parameters is θ(k−1) and an updated value of the estimates is produced as θ(k) .
By iterating the procedure K times, the resulting θ(k) is the maximizer for an approximate Algorithm 1 A generic iteration of SAEM-ABC using acceptance-rejection Simulation step: here we update X (k) using an ABC procedure sampling from π δ (X|Y; θ(k−1) ): Repeat -Generate a candidate X * from the latent model dynamics conditionally on θ(k−1) , either by numerical methods or using the transition density (if available) i.e. by generating using the exact law p X (•; θ(k−1) ). - likelihood (the approximation implied by using ABC).The "repeat loop" can be considerably expensive using a distance ρ(Y * , Y) ≤ δ as acceptance of Y * (hence acceptance of X * ) is a rare event for δ reasonably small.If informative low-dimensional statistics η(•) are available, it is recommended to consider ρ(η(Y * ), η(Y)) instead.However algorithm 1 is not appropriate for state-space models, because the entire candidate trajectory X * is simulated blindly to data (an alternative approach is considered in section 3).If we consider for a moment X * ≡ X * 0:N as a generic unknown, ideally we would like to sample from the posterior π(X * |Y), which is proportional to for a given "prior" distribution π(X * ).For some models sampling from such posterior is not trivial, for example when X is a stochastic process and in that case sequential Monte Carlo methods can be used as described in section 3. A further layer of approximation is introduced when π(X * |Y) is analytically "intractable", and specifically when f (Y|X) is unavailable in closed form (though we always assume f (•|•) to be known and that it is possible to evaluate it pointwise) or is computationally difficult to approximate but it is easy to sample from.In this case ABC methodology turns useful, and an ABC approximation to π(X * |Y) can be written as Here J δ (•) is some function that depends on δ and weights the intractable posterior based on simulated data π(X * |Y * ) with high values in regions where Y and Y * are similar; therefore we would like (i) J δ (•) to give higher rewards to proposals X * corresponding to Y * having values close to Y.In addition (ii) J δ (Y, Y * ) is assumed to be a constant when Y * = Y (i.e. when δ = 0) so that J δ is absorbed into the proportionality constant and the exact marginal posterior π(X * |Y) is recovered.Basically the use of ( 6) means to simulate X * from its prior (the product of transition densities), then plug such draw into f (•|X * ) to simulate Y * , so that X * will be weighted by J δ (Y, Y * ).A common choice for J δ (•) is the uniform kernel where ρ(Y, Y * ) is some measure of closeness between Y * and Y and I is the indicator function.
A further popular possibility is the Gaussian kernel where denotes transposition, so that J δ (Y, Y * ) gets larger when Y * ≈ Y. However one of the difficulties is that, in practice, δ has to be set as a compromise between statistical accuracy (a small positive δ) and computational feasibility (δ not too small).Notice that a proposal's acceptance can be significantly enhanced when the posterior ( 6) is conditional on summary statistics of data η(Y), rather than Y itself, and in such case we would consider ρ(η(Y), η(Y * )).However, in practice for dynamical models it is difficult to identify "informative enough" summary statistics η(•), but see Martin et al. [2016] and Picchini and Forman [2015].Another important problem with the strategy outlined above is that "blind simulation" for the generation of the entire time series X * is often poor.In fact, even when the current value of θ is close to its true value, proposed trajectories rarely follow measurements when (a) the dataset is a long time series and/or (b) the model is highly erratic, for example when latent dynamics are expressed by a stochastic differential equation (section 4.2).For these reasons, sequential Monte Carlo (SMC) [Cappe et al., 2007] methods have emerged as the most successful solution for filtering in non-linear non-Gaussian state-space models.
Below we consider the ABC-SMC methodology from Jasra et al. [2012], which proves effective for state-space models.

SAEM coupled with an ABC-SMC algorithm for filtering
Here we consider a strategy for filtering that is based on an ABC version of sequential Monte Carlo sampling, as presented in Jasra et al. [2012], with some minor modifications.The advantage of the methodology is that the generation of proposed trajectories is sequential, that is the ABC distance is evaluated "locally" for each observational time point.In fact, what is evaluated is the proximity of trajectories/particles to each data point Y j , and "bad trajectories" are killed thus preventing the propagation of unlikely states to the next observation Y j+1 and so on.For simplicity we consider the case N ≡ n, h ≡ ∆.The algorithm samples from the following target density at time and for example we could take J j,δ (y j , y * j ) = I A δ,y j (y * j ) with A δ,y j = {y * j ; ρ(η(y * j ), η(y j )) < δ} as in Jasra et al. [2012] or the Gaussian kernel (7).
The ABC-SMC procedure is set in algorithm 2 with the purpose to propagate forward M simulated states ("particles").After algorithm 2 is executed, we select a single trajectory by retrospectively looking at the genealogy of the generated particles, as explained further below.The quantity ESS is the effective sample size (e.g.Liu, 2008) often estimated as ESS({w 2 and taking values between 1 and M .When considering an indicator function for J j,δ , the ESS coincides with the number of particles having positive weight [Jasra et al., 2012].Under such choice the integer M ≤ M is a lower bound (threshold set by the experimenter) on Algorithm 2 ABC-SMC for filtering Step 0.
) and normalize weights w and go to step 1.
the number of particles with non-zero weight.However in our experiments we use a Gaussian kernel and since in the examples in section 4 we have a scalar Y j , the kernel is defined as are larger for particles having Y * (m) j ≈ Y j .We consider "stratified resampling" [Kitagawa, 1996] in step 1 of algorithm 2.
In addition to the procedure outlined in algorithm 2, once the set of weights {w n } is available at time t n , we propose to follow Andrieu et al. [2010] (see their PMMH algorithm) and sample a single index from the set {1, ..., M } having associated probabilities {w , with b m n := m .Hence, at the end of algorithm 2 we can sample m and construct its genealogy: the sequence of states {X t } resulting from the genealogy of m is the chosen path that will be passed to SAEM, see algorithm 3. The selection of this path is crucially affected by "particles impoverishment" issues, see below.[An alternative procedure, which we do not pursue here, is to sample at time t n not just a single index m , but instead sample with replacement say G ≥ 1 times from the set {1, ..., M } having associated probabilities {w n }.Then construct the genealogy for each of the G sampled indeces, and for each of the resulting G sampled paths X g,k calculate the corresponding vector-summaries S g,k c := S c (Y, X g,k ), g = 1, ..., G. Then it would be possible to take the sample average of those G summaries Sk c as in Kuhn and Lavielle [2005], and plug this average in place of S c (Y, X (k) ) in step 3 of algorithm 3. Clearly this approach increases the computational complexity linearly with G.] Notice that in Jasra et al. [2012] n ABC thresholds {δ 1 , ..., δ n } are constructed, one threshold for each corresponding sampling time in {t 1 , ..., t n }: these thresholds do not need to be set by the user but can be updated adaptively using a stochastic data-driven procedure.This is possible because the ABC-SMC algorithm in Jasra et al. [2012] is for filtering only, that is θ is a fixed known quantity.However in our scenario θ is unknown, and letting the thresholds vary adaptively (and randomly) between each pair of iterations of a parameter estimation algorithm is not appropriate.This is because the evaluation of the likelihood function at two different iterations k and k of SAEM would depend on a procedure determining corresponding (stochastic) sequences {δ n }.Therefore the likelihood maximization would be affected by the random realizations of the sequences of thresholds.In our case, we let the threshold vary deterministically: that is we choose a δ common to all time points {t 1 , ..., t n } and execute a number of SAEM iterations using such threshold.Then we deterministically decrease the threshold according to a user-defined scheme and execute further SAEM iterations, and so on.Our SAEM-ABC procedure is detailed in algorithm 3 with a user-defined sequence Here K is the number of SAEM iterations, as defined at the end of section 2.1.In our applications we show how the algorithm is not overly sensitive to small variations in the δ's.
Regarding the choice of δ values in applications, recall that the interpretation of δ in equation ( 8) is that of the standard deviation for a perturbed model.This implies that a synthetic observation at time t j , denoted with Y * j , can be interpreted as a perturbed version of Y j , where the observed Y j is assumed generated from the state-space model in equation ( 1), while Y * j ∼ N (Y j , δ 2 ).With this fact in mind, δ can easily be chosen to represent some deviation from the actual observations.Therefore, typically it is enough to look at the time evolution of the data, to guess at least the order of magnitude of the starting value for δ.
Finally, in our applications we compare SAEM-ABC with SAEM-SMC.SAEM-SMC is detailed in algorithm 4: it is structurally the same as algorithm 3 except that it uses the bootstrap filter to select the state trajectory.The bootstrap filter [Gordon et al., 1993] is just like algorithm 2, except that no simulation from the observations equation is performed (i.e. the Y * (m) j are not generated) and , hence there is no need to specify the δ.The trajectory X (k) selected by the bootstrap filter at kth SAEM iteration is then used to update the statistics s k .
As studied in detail in section 4.3, the function J t,δ has an important role in approximating the density π(X This means that for a small δ (which is set by the user) J δ assigns large weights only to very promising particles, that is those particles where N (•|a; b) is a Gaussian distribution with mean a and variance b, then the bootstrap filter underlying SAEM-SMC assigns weights proportionally to the measurements density f (Y t |X t ; σ ε ), which is affected by the currently available (and possibly poor) value of σ ε , when σ ε is one of the unknowns to be estimated.When σ ε is overestimated, as in section 4.3, there is a risk to sample particles which are not really "important".Since trajectories selected in step 2 of algorithm 4 are drawn from π(X 1:t |Y 1:t ), clearly the issues just highlighted contribute to bias the inference.
A further issue is studied in section 4.1.3.There we explain how SAEM-ABC, despite being an approximate version of SAEM-SMC (due to the additional approximation induced by using a strictly positive δ) can in practice outperform SAEM coupled with the simple bootstrap filter, because of "particles impoverishment" problems.Particles impoverishment is a pathology due to frequently implementing particles resampling: the resampling step reduces the "variety" of the particles, by duplicating the ones with larger weights and killing the others.We show that when particles fail to get close to the targeted observations, then resampling is frequently triggered, this degrading the variety of the particles.However with an ABC filter, particles receive some additional weighting, due to the function J j,δ in equation ( 8) taking values J j,δ > 1 for small δ when Y * j ≈ Y j , which allows a larger number of particles to have a non-negligible weight and therefore different particles are resampled, this increasing their variety (at least for the application in section 4.1.1,but this is not true in general).This is especially relevant in early iterations of SAEM, where θ (k) might be far from its true value and therefore many particles Algorithm 3 SAEM-ABC using a particle filter Step 0. Set parameters starting values θ(0) , set M , M and k := 1. Set the sequence {δ 1 , ..., δ L } and δ := δ 1 .
Step 1.For fixed θ(k−1) apply the ABC-SMC algorithm 2 with threshold δ, M particles and particles threshold M . 2 Sample an index m from the probability distribution {w Step 0. Set parameters starting values θ(0) , set M , M and k := 1.For fixed θ(k−1) apply the bootstrap filter below using M particles and particles threshold M .
Step 2. Sample an index m from the probability distribution {w n } on {1, ..., M } and form the path X (k) resulting from the genealogy of m .
Step 3. Stochastic Approximation step : update of the sufficient statistics Set k := k + 1. Go to step 1a.might end far from data.By letting δ decrease not "too fast" as SAEM iterations increase, we allow many particles to contribute to the states propagation.However, as δ approaches a small value, only the most promising particles will contribute to selecting the path sampled in step 2 of algorithm 3.
These aspects are discussed in greater detail in sections 4.1.3and 4.3.However, it is of course not true that an ABC filter is in general expected to perform better than a non-ABC one.In the second example, section 4.2, an adapted (not "blind to data", i.e. conditional to the next observation) particle filter is used to treat the specific case of SDE models requiring numerical integration [Golightly and Wilkinson, 2011].The adapted filter clearly outperforms both SAEM-SMC and SAEM-ABC when the number of particles is very limited (M = 30).
Notice, our SAEM-ABC strategy with a decreasing series of thresholds shares some similarity with tempering approaches.For example Herbst and Schorfheide [2017] artificially inflate the observational noise variance pertaining f (y t |X t , •), so that the particle weights have lower variance hence the resulting filter is more stable.More in detail, they construct a bootstrap filter where particles are propagated through a sequence of intermediate tempering steps, starting from an observational distribution with inflated variance, and then gradually reducing the variance to its nominal level.

Fisher Information matrix
The SAEM algorithm allows also to compute standard errors of the estimators, through the approximation of the Fisher Information matrix.This is detailed below, however notice that the algorithm itself advances between iterations without the need to compute such matrix (nor the gradient of the function to maximize).The standard errors of the parameter estimates can be calculated from the diagonal elements of the inverse of the Fisher information matrix.Its direct evaluation is difficult because it has no explicit analytic form, however an estimate of the Fisher information matrix can easily be implemented within SAEM as proposed by Delyon et al. [1999] using the Louis' missing information principle [Louis, 1982].
The Fisher information matrix (θ) = L(Y; θ) can be expressed as: where denotes transposition.An on-line estimation of the Fisher information is obtained using the stochastic approximation procedure of the SAEM algorithm as follows (see Lavielle, 2014 for an off-line approach).At the (k + 1)th iteration of the algorithm, we evaluate the three following quantities: As the sequence ( θ(k) ) k converges to the maximum of an approximate likelihood, the sequence (F k ) k converges to the corresponding approximate Fisher information matrix.It is possible to initialize G 0 and H 0 to be a vector and a matrix of zeros respectively.We stress that we do not make use of the (approximate) Fisher information during the optimization.

Simulation studies
Simulations were coded in MATLAB (except for R examples using the pomp package) and executed on a Intel Core i7-2600 CPU 3.40 GhZ.Software code can be found online either at https://github.com/umbertopicchini/SAEM-ABCor as supplementary material in the version of this paper published on Computational Statistics (doi:10.1007/s00180-017-0770-y),see also the footnote on page 2. For all examples we consider a Gaussian kernel for J j,δ as in (8).
As described at the end of section 2.1, in SAEM we always set γ k = 1 for the first K 1 iterations and γ k = (k − K 1 ) −1 for k ≥ K 1 as in Lavielle [2014].All results involving ABC are produced using algorithm 3 i.e. using trajectories selected via ABC-SMC.We compare our results with standard algorithms for Bayesian and "classical" inference, namely Gibbs sampling and particle marginal methods (PMM) [Andrieu and Roberts, 2009] for Bayesian inference and the improved iterated filtering (denoted in literature as IF2) found in Ionides et al. [2015] for maximum likelihood estimation.In order to perform a fair comparison between methods, we make use of well tested and maintained code to fit models with PMM and IF2 via the R pomp package [King et al., 2015].All the methods mentioned above use sequential Monte Carlo algorithms (SMC), and their pomp implementation considers the bootstrap filter.We remark that our goal is not to consider specialized state-of-art SMC algorithms, with the notable exception mentioned below.
Our focus is to compare the several inference methods above, while using the bootstrap filter for the trajectory proposal step: the bootstrap filter is also the approach considered in King et al. [2015] and Fasiolo et al. [2016] hence it is easier for us to compare methods using available software packages such as pomp.However in section 4.2 we also use a particles sampler that conditions upon data and is suitable for state-space models having latent process expressed by a stochastic differential equation [Golightly and Wilkinson, 2011].

Non-linear Gaussian state-space model
Consider the following Gaussian state-space model with ν j , τ j ∼ N (0, 1) i.i.d. and X 0 = 0. Parameters σ x , σ y > 0 are the only unknowns and therefore we conduct inference for θ = (σ 2 x , σ 2 y ).We first construct the set of sufficient statistics corresponding to the complete log-likelihood L c (Y, X).This is a very simple task since Y j |X j ∼ N (X j , σ 2 y ) and X j |X j−1 ∼ N (2 sin(e X j−1 ), σ 2 x ) and therefore it is easy to show that S σ 2 x = n j=1 (X j − 2 sin(e X j−1 )) 2 and S σ 2 y = n j=1 (Y j − X j ) 2 are sufficient for σ 2 x and σ 2 y respectively.By plugging these statistics into L c (Y, X) and equating to zero the gradient of L c with respect to (σ 2 x , σ 2 y ), we find that the M-step of SAEM results in updated values for σ 2 x and σ 2 y given by S σ 2 x /n and S σ 2 y /n respectively.Expressions for the first, second and mixed derivatives, useful to obtain the Fisher information as in Section 3, are given in appendix.

Results
We generate n = 50 observations for {Y j } with σ 2 x = σ 2 y = 5, see Figure 1.We first describe results obtained using SAEM-ABC.Since the parameters of interest are positive, for numerical convenience we work on the log-transformed versions (log σ x , log σ y ).Our setup consists in running 30 independent experiments with SAEM-ABC: for each experiment we simulate parameter starting values for (log σ x , log σ y ) independently generated from a bivariate Gaussian distribution with mean the true value of the parameters, i.e. (log √ 5, log √ 5), and diagonal covariance matrix having values (2,2) on its diagonal.For all 30 simulations we use the same data and the same setup except that in each simulation we use different starting values for the parameters.For each of the 30 experiments we let the threshold δ decrease in the set of values δ ∈ {2, 1.7, 1.3, 1} for a total of K = 400 SAEM-ABC iterations, where we use δ = 2 for the first 80 iterations, δ = 1.7 for further 70 iterations, δ = 1.3 for further 50 iterations and δ = 1 for the remaining 200 iterations.The influence of this choice is studied below.As explained in section 3, the largest value for δ can be set intuitively, by looking at Figure 1, where it is apparent that considering deviations δ = 2 of the simulated observations Y * j from the actual observation Y j should be reasonable.For example the empirical standard deviation of the differences |Y j − Y j−1 | is about 2. Then we let δ decrease progressively as SAEM-ABC evolves.We take K 1 = 300 as the number of SAEM warmup iterations and use different numbers of particles M in our simulation studies, see Table 1.We impose resampling when the effective sample size ESS gets smaller than M , for any attempted value of M .At first we show that taking M = 200 gives good results for SAEM-ABC but not for SAEM-SMC, see Table 1.Table 1 reports the median of the 30 estimates and their 1 st − 3 rd quartiles: we notice that M = 1, 000 particles are able to return satisfactory estimates when using SAEM-ABC.Figure 2 shows the rapid convergence of the algorithm towards the true parameter values for all the 30 repetitions (though difficult to notice visually, several simulations start at locations very far from the true parameter values).Notice that it only required about 150 seconds to perform all 30 simulations: this is the useful return out of the effort of constructing the analytic quantities necessary to run SAEM.Also, the algorithm is not very sensitive to the choice of δ s.For example, we also experimented with δ ∈ {4, 3, 2, 1} and we obtained very similar results, for example our thirty experiments with (M, M ) = (1000, 200) resulted in medians (1st-3rd quartiles) σx = 2.30 [2.27,2.35],σy = 1.90 [1.86,1.91],which are very close to the ones in Table 1.Finally, notice that results are not overly sensitive to the way δ is decreased: for example, if we let δ decrease uniformly with steps of size 1/3, that is δ ∈ {2, 1.67, 1.33, 1} with δ = 2 for the first 50 iterations, then let it decrease every 50 iterations until δ = 1, we obtain (when M = 1, 000 and M = 200) σx = 2.30 [2.27,2.33],σy = 1.91 [1.87,1.95],compare with Table 1.
We then perform 30 simulations with SAEM-SMC using the same simulated data and parameters starting values as for SAEM-ABC.As from Table 1 simulations for σ y converge to completely wrong values.For this case we also experimented with M = 5, 000 using M = 2, 000 but this does not solve the problem with SAEM-SMC, even if we let the algorithm start at the true parameter values.We noted that when using M = 2, 000 particles with SAEM-ABC and M = 200 the algorithm resamples every fourth observation, and in a generic iteration we observed an ESS of about 100 at the last time point.Under the same setup SAEM-SMC resampled at each time point and resulted in an ESS of about 10 at the last time point (a study on the implications of frequent resampling is considered in section 4.1.3).Therefore we now perform a further simulation study to verify whether using a smaller M (hence reducing the number of times resampling is performed) can improve the performance of SAEM-SMC.Indeed, using for example M = 20 gives better results for SAEM-SMC, see Table 1 (there we only use M = 1, 000 for comparison between methods).This is further investigated in section 4.1.3.

Comparison with iterative filtering and a pseudo-marginal Bayesian algorithm
We compare the results above with the iterated filtering methodology for maximum likelihood estimation (IF2, Ionides et al., 2015), using the R package pomp [King et al., 2015].We provide pomp with the same data and starting parameter values as considered in SAEM-ABC and SAEM-SMC.We do not provide a detailed description of IF2 here: it suffices to say that in IF2 particles are generated for both parameters θ (e.g. via perturbations using random walks) and for the systems state (using the bootstrap filter).Moreover, same as with ABC methods, a "temperature" parameter (to use an analogy with the simulated annealing optimization method) is let decrease until the algorithm "freezes" around an approximated MLE.This temperature parameter, here denoted with , is decreased in ∈ {0.9, 0.7, 0.4, 0.3, 0.2} which seems appropriate as explained below, where the first value is used for the first 500 iterations of IF2, then each of the remaining values is used for 100 iterations, for a total of 900 iterations.Notice that the tested version of pomp (v.1.4.1.1)uses a bootstrap filter that resamples at every time point, hence results obtained with IF2 are not directly comparable with SAEM-ABC and SAEM-SMC.Results are in Table 1 and a sample output from one of the simulations obtained with M = 1, 000 is in Figure 3. From the loglikelihood in Figure 3 we notice that the last major improvement in likelihood maximization takes place at iteration 600 when becomes = 0.7.Reducing further does not give any additional benefit (we have verified this in a number of experiments with this model).
Notice that in order to run, say, 400 iterations of IF2 with M = 1, 000 for a single experiment, instead of thirty, it required about 70 seconds.That is IF2 is about fourteen times slower than SAEM-ABC, although the comparison is not completely objective as we coded SAEM-ABC with Matlab, while IF2 is provided in an R package with forward model simulation implemented in C. Finally we use a particle marginal method (PMM, Andrieu and Roberts, 2009) on a single simulation (instead of thirty), as this is a fully Bayesian methodology and results are not directly comparable with SAEM nor IF2.The PMM we construct approximates the likelihood function of the state space model using a bootstrap filter with M particles, then plugs such likelihood approximation into a Metropolis-Hastings algorithm.As remarkably shown in Beaumont [2003] and Andrieu and Roberts [2009], a PMM returns a Markov chain for θ having the posterior π(θ|Y 1:n ) as its stationary distribution.This implies that PMM is an algorithm producing exact Bayesian inference for θ, regardless the number of particles used to approximate the likelihood.
Once more we make use of facilities provided in pomp to run PMM.We set uniform priors U (0.1, 15) for both σ x and σ y and run 4,000 MCMC iterations of PMM using 2,000 particles for the likelihood approximation.Also, we set the PMM algorithm in the most favourable way, by starting it at the true parameter values (we are only interested in the inference results, rather than showing the performance of the algorithm when starting from remote values).The proposal function for the parameters uses an adaptive MCMC algorithm based on Gaussian random walks, and was tuned to achieve the optimal 7% acceptance rate [Sherlock et al., 2015].For this single simulation, we obtained the following posterior means and 95% intervals: σx = 1.52 [0.42,2.56],σy = 1.53 [0.36,2.34].

The particles impoverishment problem
Figure 4 reports the effective sample size ESS as a function of time t (at a generic iteration of SAEM-SMC, the 20th in this case).As expected, for a smaller value of M the ESS is most of times smaller than when a larger M is chosen, with the exception of a few peaks.This is a direct consequence of performing resampling more frequently when M is larger.However, a phenomenon that is known to be strictly linked to resampling is that of "samples impoverishment", that is the resampling step reduces the "variety" of particles, by duplicating the ones with larger weight and killing the others.In fact, when many particles have a common "parent" at time t, these are likely to end close to each other at time t + 1.This has a negative impact on the inference because the purpose of the particles is to approximate the density (5) (or ( 6)) which generates the trajectory sampled at step 2 in algorithms 3 and 4. Lack of variety in the particles reduces the quality of this approximation.Indeed, Figure 5 shows that the variety of the particles (as measured by the number of distinct particles) gets impoverished for a larger M , notice for example that the solid line in Figure 5 almost always reaches its maximum attainable value M = 1, 000, that is all particles are distinct, while this is often not the case for the dashed line.Since the trajectory X (k) that is selected at iteration k of SAEM, either in step 2 of algorithm 3 or in step 2 of algorithm 4, follows from backward-tracing the genealogy of a certain particle, having variety in the cloud of particles is crucial here.
This seems related to the counter-intuitive good performance of the ABC-filter, even though SAEM-ABC is based on the additional approximation induced by the J function in equation ( 8).We now produce plots for the ESS values and the number of distinct particles at the smallest value of the ABC threshold δ.As we see in Figure 6, while the ESS are not much different from the ones in Figure 4, instead the number of distinct particles in Figure 7 is definitely higher than the SAEM-SMC counterpart in Figure 5, meaning that such number drops below the maximum M = 1, 000 fewer times.For example from Figure 5 we can see that when ( M , M ) = (200, 1000) the number of distinct particles drops 19 times away from the maximum M = 1, 000, when using SAEM-SMC.With SAEM-ABC this number drops only 13 times (Figure 7).When the number of resampling steps is reduced, using ( M , M ) = (20, 1000), we have more even results, with the number of distinct particles dropping six times for SAEM-SMC and five times for SAEM-ABC.
As a support to this remark, refer to Table 2: we perform thirty independent estimation procedures, and in each we obtain the sample mean of the ESS (means computed over varying time t) and the sample mean of the number of distinct particles (again over varying t).Then, we report the mean over the thirty estimated sample means of the ESS (and corresponding standard deviation) and the same for the number of distinct particles.Clearly numbers are favourable to SAEM-ABC, showing consistently larger values (with small variation between the thirty experiments).
We argue that when many particles, as generated in the bootstrap filter, fail to get close to the target observations, then resampling is frequently triggered, this degrading the variety of the particles.However with an ABC filter, particles receive some additional weighting (due to a J j,δ > 1 in (8) when Y * (m) j ≈ Y j and a small δ) which allows for a larger number of particles to have a non-negligible weight.While it is not true that an ABC filter is in general expected to perform better than a non-ABC one, here we find that the naive bootstrap filter performs worse than the ABC counterpart for the reasons discussed above.

Relation with Gibbs sampling
As previously mentioned, the generation of s k in (3) followed by corresponding parameter estimates θ(k) in ( 4) is akin to two steps of a Gibbs sampling algorithm, with the important distinction that in SAEM θ(k) is produced by a deterministic step (for given s k ).We show that the construction of a Gibbs-within-Metropolis sampler is possible, but that a naive implementation fails while a "non-central parametrization" seems necessary [Papaspiliopoulos et al., 2007].For given initial vectors (X (0) , σ (0) x , σ (0) y ), we alternate sampling from the conditional distributions (i)   p(X (b) |σ x , X (b) , Y) where b represents the iteration counter in the Gibbs sampler (b = 1, 2, ...).The resulting multivariate sample (X (b) , σ y ) is a draw from the posterior distribution π(X, σ x , σ y |Y).We cannot sample from the conditional densities in (i)-(iii), however at a generic iteration b it is possible to incorporate a single Metropolis-Hastings step targeting the corresponding densities in (i)-(iii) separately, resulting in a Metropolis-within-Gibbs sampler, see e.g.Liu [2008].Notice that what (i) implies is a joint sampling (block-update) for all the in X, however it is also possible to sample individual elements one at a time (single-site update) from p(X For both single-site and block-update sampling, the mixing of the resulting chain is very poor, this resulting from the high correlations between sampled quantities, notably the correlation between elements in X and also between X and (σ x , σ y ).However, there exists a simple solution based on breaking down the dependence between some of the involved quantities: for example at iteration b we propose in block a vector X # generated "blindly" (i.e.conditionally on the current parameter values (σ ), but unconditionally on Y) by iterating through (9) then accept/reject the proposal according to a Metropolis-Hastings step, then define X * (b) := X/σ , where X is the last accepted proposal of X (that is X ≡ X # if X # has been accepted).Sample from (ii) p(σ x , X * (b) , Y), and transform back to X := σ (b) x X * (b) .This type of updating scheme is known as non-central parametrization [Papaspiliopoulos et al., 2007].The expressions for the conditional densities in (i)-(iii) are given in the appendix for the interested reader.Note that the SAEM-ABC and SAEM-SMC do not require a non-central parametrization.The maximization step of SAEM smooths out the correlation between the proposed parameter and the latent state, and the numerical convergence of the algorithms still occur (this has also been noticed in previous papers on SAEM-MCMC, see Kuhn andLavielle, 2005, Donnet andSamson, 2008).
Once more, we attempt estimating data produced with (σ x , σ y ) = (2.23,2.23) (as in section 4.1.1)this time using Metropolis-within-Gibbs. Trace plots in Figures 8 show satisfactory mixing for three chains starting at three random values around (σ x , σ y ) = (6, 6).However, while the chains initialized at values far from the truth rapidly approach the true values, the 95% posterior intervals fail to include them, see Figure 9.If we re-execute the same experiment, this time with data generated with smaller noise (σ x , σ y ) = (0.71, 0.71), it seems impossible to catch the ground truth parameters.From a typical chain, we obtain the following posterior means and 95% posterior intervals: σx = 1.59 [1.18,2.02],σy = 1.63 [1.28,2.08].Clearly for noisy time-dependent data, such as state-space models, particles based methods seem to better address cases where noisy data are affected by both measurement and systemic noise.

A pharmacokinetics model
Here we consider a model for pharmacokinetics dynamics.For example we could formulate a model to study the Theophylline drug pharmacokinetics.This example has often been described in literature devoted to longitudinal data modelling with random parameters (mixed-effects models), see Pinheiro and Bates [1995] and Donnet and Samson [2008].Same as in Picchini [2014] here we do not consider a mixed-effects model.We denote with X t the level of drug concentration in blood at time t (hrs).Consider the following non-authonomous stochastic differential equation: where Dose is the known drug oral dose received by a subject, K e is the elimination rate constant, K a the absorption rate constant, Cl the clearance of the drug and σ the intensity of the intrinsic stochastic noise.We simulate data at n = 100 equispaced sampling times {t 1 , t ∆ , ..., t 100∆ } = {1, 2, ..., 100} where ∆ = t j −t j−1 = 1.The drug oral dose chosen to be 4 mg.After the drug is administered, we consider as t 0 = 0 the time when the concentration first reaches X t 0 = X 0 = 8.The error model is assumed to be linear, Y j = X j + ε j where the ε j ∼ N (0, σ 2 ε ) are i.i.d., j = 1, ..., n.Inference is based on data {Y 1 , ..., Y n } collected at corresponding sampling times.Parameter K a is assumed known, hence parameters of interest are θ = (K e , Cl, σ 2 , σ 2 ε ) as X 0 is also assumed known.
Equation ( 10) has no available closed-form solution, hence simulated data are created in the following way.We first simulate numerically a solution to (10) using the Euler-Maruyama discretization with stepsize h = 0.05 on the time interval [t 0 , 100].The Euler-Maruyama scheme is given by where the {Z t } are i.i.d.N (0, 1) distributed.The grid of generated values X 0:N is then linearly interpolated at sampling times {t 1 , ..., t 100 } to give X 1:n .Finally residual error is added to X 1:n according to the model Y j = X j + ε j as explained above.Since the errors ε j are independent, data {Y j } are conditionally independent given the latent process {X t }.
Sufficient statistics for SAEM The complete likelihood is given by p(Y, X 0:N ; θ) = p(Y|X 0:N ; θ)p(X 0:N where the unconditional density p(x 0 ) is disregarded in the last product since we assume X 0 deterministic.Hence the complete-data loglikelihood is Here p(y j |x j ; θ) is a Gaussian with mean x j and variance σ 2 ε .The transition density p(x i |x i−1 ; θ) is not known for this problem, hence we approximate it with the Gaussian density induced by the Euler-Maruyama scheme, that is Notice the Gaussian distribution implied by ( 11) shares some connection with tools developed for optimal states predictions in signal processing, such as the unscented Kalman filter, e.g.Sitz et al. [2002].We now derive sufficient summary statistics for the parameters of interest, based on the complete loglikelihood.Regarding σ 2 ε this is trivial as we only have to consider n j=1 log p(y j |x j ; θ) to find that a sufficient statistic is S σ 2 ε = n j=1 (y j − x j ) 2 .Regarding the remaining parameters we have to consider N i=1 log p(x i |x i−1 ; θ).For σ 2 a sufficient statistic is Regarding K e and Cl reasoning is a bit more involved: we can write The last equality suggests a linear regression approach E(V ) = β 1 C 1 + β 2 C 2 for "responses" V i = (x i − x i−1 )/ √ x i−1 and "covariates" and β 1 = K e /Cl, β 2 = K e .By considering the design matrix C with columns C 1 and C 2 , that is C = [C 1 , C 2 ], from standard regression theory we have that β = (C C) −1 C V is a sufficient statistic for β = (β 1 , β 2 ), where denotes transposition.We take S Ke := β2 also to be used as the updated value of K e in the maximisations step of SAEM.Then we have that β1 is sufficient for the ratio K e /Cl and use β2 / β1 as the update of Cl in the M-step of SAEM.The updated values of σ and σ ε are given by S σ 2 /N and S σ 2 ε /n respectively.

Results
We consider an experiment where 50 datasets of length n = 100 each are independently simulated using parameter values (K e , K a , Cl, σ, σ ε ) = (0.05, 1.492, 0.04, 0.1, 0.1).All results pertaining SAEM-ABC use the Gaussian kernel (8).For SAEM-ABC we use a number of "schedules" to decrease the threshold δ.One of our attempts decreases the threshold as δ ∈ {0.5, 0.2, 0.1, 0.03} (results for this schedule are reported as SAEM-ABC(1) in Table 3): same as in the previous application, all we need is to determine an appropriate order of magnitude for the largest δ.
As an heuristic, we have that the empirical standard deviation of the differences |Y j − Y j−1 | is about 0.3.The first value of δ is used for the first 80 iterations then it is progressively decreased every 50 iterations.For both SAEM-ABC and SAEM-SMC we use K 1 = 250 and K = 300 and optimization started at parameter values very far from the true values: starting values are K e = 0.80, Cl = 10, σ = 0.14 and σ ε = 1.We first show results using (M, M ) = (200, 10).We start with SAEM-ABC, see Figure 10 where the effect of using a decreasing δ is evident, especially on the trajectories for σ ε .All trajectories but a single erratic one converge towards the true parameter values.Estimation results are in Table 3, giving results for three different decreasing schedules for δ, reported as SAEM-ABC(0) for δ ∈ {0.5, 0.2, 0.1, 0.05, 0.01}, the already mentioned SAEM-ABC(1) having δ ∈ {0.5, 0.2, 0.1, 0.03} and finally SAEM-ABC(2) where δ ∈ {1, 0.4, 0.1}.Results are overall satisfactory for all parameters but σ, which remains unidentified for all attempted SAEM algorithms, including those discussed later.The benefit of decreasing δ to values close to zero are noticeable, that is among the algorithms using ABC, SAEM-ABC(0) gives the best results.
From results in Table 3 and Figure 11, again considering the case M = 200, we notice that SAEM-SMC struggles in identifying most parameters with good precision.As discussed later on, when showing results with M = 30, it is clear that when a very limited amount of particles is available (which is of interest when it is computationally demanding to forward-simulate from a complex model) SAEM-SMC is suboptimal compared to SAEM-ABC at least for the considered example.However results improve noticeably for SAEM-SMC as soon as the number of particles is enlarged, say to M = 1, 000.In fact for M = 1, 000 the inference results for SAEM-ABC(0) and SAEM-SMC are basically the same.It is of course interesting to uncover the reason why SAEM-ABC(0) performs better than SEM-SMC when the number of particles is small, e.g.M = 200 (see later on for results with M = 30).In Figure 12 we compare approximations to the distribution π(X t |Y 1:t−1 ), that is the distribution of the state at time t given previous data, and the filtering distribution π(X t |Y 1:t ) including the most recent data.The former one, π(X t |Y 1:t−1 ), is approximated by kernel smoothing applied on the particles X (m) t .The filtering distribution π(X t |Y 1:t ) is approximated by kernel smoothing assigning weights W (m) t to the corresponding particles.It is from the (particles induced) discrete distribution approximating π(X t |Y 1:t ) that particles are sampled when propagating to the next time point.Both densities are computed from particles obtained at the last SAEM iteration, K, for different values of t, at the beginning of the observational interval (t = 10), at t = 40 and towards the end of the observational interval (t = 70).Since in Figure 12 we also report the true values of X t , we can clearly see that, for SAEM-SMC, as t increases the true value of X t is unlikely under π(X t |Y 1:t ) (the orange curve).In particular, by looking at the orange curve in panel (f) in Figure 12 we notice that many of the particles ending-up far from the true X t = 1.8 receive non-negligible weight.Instead in panel (e) only particles very close to the true value of the state receive considerable weight (notice the different scales on the abscissas for panel (e) and (f)).Therefore it is not unlikely that for SAEM-SMC several "remote" particles are resampled and propagated.The ones resampled in SAEM-ABC received a weight W (m) t ∝ J t,δ which is large only if the simulated observation is very close to the actual observation, since δ is very small.Therefore for SAEM-SMC, when M is small, the path sampled in step 2 of algorithm 4 (which is from π(X 1:tn |Y 1:tn )) is poor and the resulting inference biased.For a larger M (e.g.M = 1, 000) SAEM-SMC enjoys a larger number of opportunities for particles to fall close to the targeted observation hence an improved inference.
An improvement over the bootstrap filter used in SAEM-SMC is given by a methodology where particles are not proposed from the transition density of the latent process (the latter proposes particles "blindly" with respect to the next data point).For example, Golightly and Wilkinson [2011] consider a proposal distribution based on a diffusion bridge, conditionally to observed data.Their methodology is specific for state-space models driven by a stochastic differential equation whose approximate solution is obtained via the Euler-Maruyama discretisation, which is what we require.We use their approach to "propagate forward" particles.We write SAEM-GW to denote a SAEM algorithm using the proposal sampler in Golightly and Wilkinson [2011].By looking at Table 3 we notice the improvement over the simpler SAEM-SMC when M = 200.In fact SAEM-GW gives the best results of all SAEM-based algorithms we attempted,  (M, M ) = ( 200  the downside being that such sampler is not very general and is only applicable to a specific class of models, namely (i) state-space models having an SDE that has to be numerically solved via Euler-Maruayama, (ii) observations having additive Gaussian noise, and (iii) observations having a state entering linearly, e.g.y = a • x + ε for some constant a. Same as in section 4.1, we consider Bayesian estimation using a particle marginal method (PMM).PMM is run with 1000 particles for 2,000 MCMC iterations.However it turns out that PMM cannot be initialized at the same remote starting values we used for SAEM, as the approximated log-likelihood function at the starting parameter results not-finite for the considered number of particles.Therefore we let PMM start at K e = 0.05, Cl = 0.04, σ = 0.2, σ ε = 0.3 and use priors K e ∼ U (0.01, 0.2), Cl ∼ U (0.01, 0.2), σ ∼ U (0.01, 0.3), σ ε ∼ U (0.05, 0.5).
We now explore whether the ABC approach can be of aid when saving computational time is essential, for example when simulating from the model is expensive.Although the model here considered can be simulated relatively quickly (it requires numerical integration hence computing times are affected by the size of the integration stepsize h), assuming this is not the case we explore what could happen if we can only afford running a simulation with M = 30 particles.We run 100 simulations with this setting.To ease graphic representation in the presence of outliers, estimates are reported on log-scales in the boxplots in Figures 13-15.We notice that, with such a small number of particles, SAEM-ABC is still able to estimate K e and Cl accurately whereas SAEM-SMC returns more biased estimates for all parameters.Also, while both methodologies fail in estimating σ and σ ε when M = 30 SAEM-ABC is still better than SAEM-SMC confirming the previous finding that, should the computational budget be very limited, SAEM-ABC is a viable option.SAEM-GW is instead able to estimate also the residual error variability σ ε .

Summary
We have introduced a methodology for approximate maximum likelihood estimation of the parameters in state-space models, incorporating an approximate Bayesian computation (ABC) strategy.The general framework is the stochastic approximation EM algorithm (SAEM) of Delyon et al. [1999], and we embed a sequential Monte Carlo ABC filter into SAEM.SAEM requires model-specific analytic computations, at the very least the derivation of sufficient statistics for the complete log-likelihood, to approach the parameters maximum likelihood estimate with minimal computational effort.However at any iteration of SAEM, it is required the availability of a filtered trajectory of the latent systems state, which we provide via the ABC filter of Jasra et al. [2012].We call this algorithm SAEM-ABC.An advantage of using the ABC filter is its flexibility, as it is possible to modify its setup to influence the weighting of the particles.In other words, it is possible to tune a positive tolerance δ to enhance the importance of those particles that are the closest to the observations.This produced sampled trajectories for step 2 of algorithm 3 that resulted in a less biased parameter inference.This observation turned especially true for experiments run with a limited number of particles (M = 30 or 200), which is relevant for computationaly intensive models not allowing for the propagation of a large number of particles.We compared our SAEM-ABC algorithm with a version of SAEM employing the bootstrap filter of Gordon et al. [1993].The bootstrap filter is the simplest sequential Monte Carlo algorithm, and is typically a default option in many software packages (e.g. the smfsb and pomp R packages, Wilkinson, 2015 andKing et al., 2015 respectively).In our work we show that, in some cases, SAEM-SMC (that is SAEM using a bootstrap filter) is sometimes inferior to SAEM-ABC, for example when the number of particles is small (section 4.3) or when too frequent resampling causes "particles impoverishment" (section 4.1.3).SAEM-ABC requires the user to specify a sequence of ABC thresholds δ 1 > • • • δ L > 0. For one-dimensional time-series setting these thresholds is intuitive, since these represent standard deviations of perturbed (simulated) observations.Therefore the size of the largest one (δ 1 ) can be determined by looking at plots of the observed time-series.
Ultimately, while we are not claiming that an ABC filter should in general be preferred to a non-ABC filter, as the former one induces some approximation, it can be employed when it is difficult to construct (or implement) a more advanced sequential Monte Carlo filter.In our second application we consider a sequential Monte Carlo filter due to Golightly and Wilkinson [2011].This filter (which we call GW) is specifically designed for state-space models driven by stochastic differential equations (SDEs) requiring numerical discretization.While GW improves noticeably over the basic bootstrap filter, GW is not very general: again, it is specific for statespace models driven by SDEs; the observation equation must have latent states entering linearly and measurement errors must be Gaussian distributed.The ABC approach instead does not impose any limitation on the model structure.
The four coordinates of the gradient are: Entries for the Fisher information matrix are (recall this is a symmetric matrix, therefore redundant terms are not reported.Further missing entries consist of zeros): Denote with m such index and with a m j the "ancestor" of the generic mth particle sampled at time t j+1 , with 1 ≤ a m j ≤ M (m = 1, ..., M , j = 1, ..., n).Then we have that particle m has ancestor a m n−1 and in general particle m at time t j+1 has ancestor b m j n } on {1, ..., M } and form the path X(k) resulting from the genealogy of m .Step 3. Stochastic Approximation step : update of the sufficient statistics

Figure 8 :
Figure 8: Nonlinear Gaussian model: three chains with different starting values from the Metropolis-within-Gibbs sampler for data generated with (σx, σy) = (2.23,2.23).Horizontal lines are the true parameter values.

Figure 13 :Figure 14 :Figure 15 :
Figure 13: Theophylline model: estimates obtained with SAEM-SMC when M = 30.From left to right: log Ke, log Cl, log σ and log σε.Green lines are the true parameter values on log-scale.

Table 1 :
Non-linear Gaussian model: medians and 1 st − 3 rd quartiles for estimates obtained on 30 independent simulations.(*)The IF2 method resamples at every time point, while SAEM-ABC and SAEM-SMC resamples only when ESS < M .

Table 2 :
Nonlinear -Gaussian model: mean ESS and mean number of distinct particles (and corresponding standard deviations in parentheses) at a generic iteration of SAEM-ABC and SAEM-SMC.Averages and standard deviations are taken over 30 independent repetitions of the experiment.

Table 3 :
Theophylline: medians and 1 st − 3 rd quartiles for estimates obtained on 50 independent simulations