Pierre-Arnaud Coquelin

Romain Deguest

Rémi Munos

Numerical methods for sensitivity analysis of Feynman-Kac models

The aim of this work is to provide efficient numerical methods to estimate the gradient of a Feynman-Kac flow with respect to a parameter of the model. The underlying idea is to view a Feynman-Kac flow as an expectation of a product of potential functions along a canonical Markov chain, and to use usual techniques of gradient estimation in Markov chains. Combining this idea with the use of interacting particle methods enables us to obtain two new algorithms that provide tight estimations of the sensitivity of a Feynman-Kac flow. Each algorithm has a linear computational complexity in the number of particles and is demonstrated to be asymptotically consistent. We also carefully analyze the differences between these new algorithms and existing ones. We provide numerical experiments to assess the practical efficiency of the proposed methods and explain how to use them to solve a parameter estimation problem in Hidden Markov Models. To conclude we can say that these algorithms outperform the existing ones in terms of trade-off between computational complexity and estimation quality.

Introduction

Let us consider a Hidden Markov Model (HMM) parameterized by θ ∈ Θ, where Θ ⊂ R n θ is an open subset of R n θ . The state process (X p,θ ) p≥0 is an homogeneous Markov chain with initial probability measure µ θ (dx0) and Markov transition kernel K θ (dxp+1|xp). The observation process (Y p,θ ) p≥1 is linked with the state process by the conditional probability measure P(Yp ∈ dyp|Xp = xp) = g θ (yp, xp)λ(dyp), where λ is a fixed probability measure. Suppose we are given a sequence of successive realizations of the observation process (y p,θ * ) p≥1 which were obtained using an unknown parameter denoted by θ * . Our objective is to recover θ * using the sequence of observations (y p,θ * ) p≥1 . In the following, we make a slight abuse of notation by writing g p,θ (xp) = g θ (y p,θ * , xp) and we adopt the usual notation zi:j = (z k ) i≤k≤j .

It is known that, under some technical conditions on HMMs [START_REF] Douc | Asymptotics of the maximum likelihood estimator for general hidden markov models[END_REF], the maximum likelihood estimator (MLE) provides an asymptotically consistent estimation of θ * . Indeed, by defining the log-likelihood function at the time step n as ln(θ) = log P(Y 1:n,θ ∈ dy 1:n,θ * ), the sequence of estimated parameters θn = argmax θ∈Θ ln(θ) converges to the true parameter θ * . Thus, the problem of parameter estimation is transposed into the problem of maximizing the log-likelihood function. Except for some particular cases such as finite state space models or linear gaussian models, ln cannot be analytically computed and requires to be numerically estimated (for example using particle filters [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF]). In these situations, maximizing the log-likelihood function may be considered as a stochastic optimization problem. A natural approach to solve it is to use a stochastic gradient ascent on ln(θ). This consists in applying a Robbin-Monroe procedure to search for a zero of the differential of the log-likelihood function ∇ln with respect to the parameter θ. The crucial point for the practical success of such an optimization procedure is to provide a tight estimation ∇ln ≈ ∇ln of the gradient of the log-likelihood function.

The initial problem of parameter estimation in HMMs may be undergone by computing a tight estimator of the gradient of the log-likelihood function. Using elementary bayesian calculus (see [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF] for details), one can show that:

ln(θ) = n p=1 log[π p|p-1,θ (g p,θ )], (1)
where we define the predicted filtering distribution π p|p-1,θ (dxp) = P(X p,θ ∈ dxp|Y 1:p-1,θ = y 1:p-1,θ * ) for p ≥ 1, and where π(g) refers to g(x)π(dx) (all precise definitions will be detailed in the Notations Section). Differentiating Equation (1) leads to:

∇ln(θ) = n p=1 ∇[π p|p-1,θ (g p,θ )] π p|p-1,θ (g p,θ ) . (2)
In order to evaluate the gradient of the log-likelihood, we need to compute the flows π p|p-1,θ (g p,θ ) and the gradient of these flows ∇[π p|p-1,θ (g p,θ )].

The flows π p|p-1,θ (g p,θ ) can be evaluated using standard algorithms as particle filters. Nonetheless, the problem of efficiently estimating the gradient of the flows ∇[π p|p-1,θ (g p,θ )] is still an open question. It is the subject of the present contribution.

Before going further, let us introduce some new notations. It is well known that π p|p-1,θ satisfies the following bayesian recurrence relation (see for example [START_REF] Doucet | On sequential monte carlo sampling methods for bayesian filtering[END_REF]):

π p+1|p,θ (dxp+1) = Xp K θ (dxp+1|xp)g θ (xp)π p|p-1,θ (dxp) Xp g θ (xp)π p|p-1,θ (dxp)
.

(

) 3
Using the Boltzman-Gibbs measure Ψg,π(dx) = g(x)π(dx) π(g)

, the relation (3) becomes π p+1|p,θ = Ψg p,θ ,π p|p-1,θ K θ .

More generally, given a sequence of positive measurable functions (Gp) p≥0 and a Markov chain with initial probability distribution µ and transition kernels (Mp) p≥1 , the sequence of marginal predicted Feynman-Kac measures (ηp) p≥0 is recursively defined by the following measure-valued equations:

η0 = µ ηp+1 = ΨG p,ηp Mp+1. (5)
Equation [START_REF] Chan | Optimisation of particle filters using simultaneous perturbation stochastic approximation[END_REF] shows that the predicted filtering distribution π p+1|p,θ is equal to the predicted Feynman-Kac measure ηp associated to the Markov chain with initial measure µ and transition kernels Mp = K θ , and to the potential functions Gp = gp.

Similarly to what has been described earlier, in the Feynman-Kac framework, the initial parameter estimation problem may be turned into the problem of estimating the sensitivity ∇ηn(f ) of a Feynman-Kac flow with respect to a parameter of a Feynman-Kac model. The framework of Feynman-Kac formulae [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF] is a subject of high interest which appears in many applications besides Hidden Markov Models, such as directed polymer simulations, random excursion models or genetic optimization algorithms. In the remainder of this paper, we adopt the general Feynman-Kac formalism. It provides more general results, a better physical intuition of measure representations via interacting particle interpretation, and it allows to use remarkable convergence results.

To the best of our knowledge, the problem of gradient estimation of flows was never addressed directly in the Feynman-Kac framework, but it is the object of a lot of efforts in the field of HMMs. Previous works can be cast into three categories: joint particle estimation of the linear tangent filter and the filter by differentiating the update and prediction equations of the filter [START_REF] Fichoud | Particle-based methods for parameter estimation and tracking : numerical experiments[END_REF][START_REF] Guyader | A particle implementation of the recursive mle for partially observed diffusions[END_REF][START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF], joint particle estimation of the linear tangent filter and the filter by differentiating the prediction equation and using a likelihood ratio approach to deal with the update step [START_REF] Doucet | Parameter estimation in general statespace models using particle methods[END_REF][START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF], and direct sensitivity analysis of the Markov chain followed by the particles [START_REF] Poyiadjis | Particle filter as a controlled markov chain for on-line parameter estimation in general state space models[END_REF][START_REF] Poyiadjis | Particle Method for Parameter Estimation in General State Space Models[END_REF]. More precisely, [START_REF] Fichoud | Particle-based methods for parameter estimation and tracking : numerical experiments[END_REF][START_REF] Guyader | A particle implementation of the recursive mle for partially observed diffusions[END_REF][START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF] address the problem of jointly representing by a particle population the linear tangent measure of the filter and the measure of the filter. A representation of the linear tangent transition kernel from which one can sample is assumed to be given. It is used to compute some linear tangent update and prediction equations. The algorithm transports a joint particle approximation of the filter and its gradient using these linear tangent equations. An extensive compilation of experimental results on linear gaussian models is available in [START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF]. In [START_REF] Doucet | Parameter estimation in general statespace models using particle methods[END_REF] the authors address the same problem as above. The filter and its gradient are represented by the same particle population, but using different weights. As the population is selected using the weights associated to the filter, the weights associated to the linear tangent filter are multiplied by likelihood ratios during selection procedures. This operation prevents the algorithm from being efficient because of the resulting additional variance term which explodes with the number of time-steps. In [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] the authors start with the same idea as in [START_REF] Doucet | Parameter estimation in general statespace models using particle methods[END_REF], but instead, they used a marginal particle filter to compute the weights associated to the particle approximation of the linear tangent filter. This approach solves the stability problem encountered in previous works. Unfortunately, the computational complexity is quadratic in the number of particles, thus preventing this method to be applied to real world applications. In [START_REF] Poyiadjis | Particle filter as a controlled markov chain for on-line parameter estimation in general state space models[END_REF][START_REF] Poyiadjis | Particle Method for Parameter Estimation in General State Space Models[END_REF], the authors proposed to approximate the sensitivity of the filter by the sensitivity of the Markov chain followed by its particles representation. Usual gradient evaluation techniques issued from controlled Markov chain theory are used. The obtained algorithms suffer from high variance and are not not of practical interest. As far as we know, at this time, [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] is the state of art method.

The aim of this work is to design efficient numerical methods for estimating the sensitivity of a Feynman-Kac flow with respect to a parameter of a the Feynman-Kac model. Our contribution consists in providing two sensitivity estimates, proven to be consistent, whose approximation quality is comparable to [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF], but whose computational complexity is only linear in the number of particles.

The paper is divided into three parts. In a first section we give a formal introduction to Feynman-Kac models in order to precisely set the gradient estimation problem. The second part is the heart of this contribution. The basic idea is to view a Feynman-Kac flow as an expectation of the product of the potential functions Gp along the canonical chain (Xp) p≥0 . Using basic calculus on this expectation provides a formal expression of the gradient of a F-K flow in the form of an integral over the F-K path. Any interacting particle method can be used to estimate tightly this integral and in the same time, to provide a tight estimation of the sensitivity of the F-K flow. The proposed algorithms have a linear computational complexity in the number of particles and a Law of Large Numbers is stated. In the third part we provide numerical experiments, on two different gradient evaluation problems in HMMs, to assess the practical efficiency of the proposed algorithms. In particular, we compare the bias, variance and CPU time of these algorithms with existing methods. We experimentally show that the approximation error does not increase with the number of time steps and that this error tends to zero when the number of particles tends to infinity. We also show, since this was our initial motivation for this work, that these algorithms can be used to solve the parameter estimation problem in nonlinear HMMs.

Notations

All random variables are defined on a common probability space (Ω, F , P). Most of random variables, measures and functions depend on a parameter θ ∈ Θ where Θ is an open subset of R n θ . Nonetheless, when there is no ambiguity, we omit to note explicitly this dependancy. A state space (E, σ(E)) is said to be general if its σ-field σ(E) is countably generated. In the following all state spaces are assumed to be general and for any d > 0, R d is equipped with the Borel σ-algebra σ(R d ). We will use the following notations:

• The gradient operator ∇ denotes the derivative with respect to θ, and the symbol denotes the derivative with respect to the state variable.

• To each state space E, we associate the following sets: D(E), the set of probability measures on E, B(E) the space of measurable functions from E to R and

B + (E) = {f ∈ B(E)|f ≥ 0, f = 0}
the space of positive measurable functions from E to R that are not everywhere null. All along the paper we work with test functions defined from E to R, but the results obtained remain the same if one works with functions defined from E to R m where m ≥ 1.

• Given two state spaces E and F , a Markov transition kernel M from E to F is an application

M : E × σ(F ) → [0, 1] such that ∀A ∈ σ(F ), {x → M (x, A)} ∈ B(E) and ∀x ∈ E, {A → M (x, A)} ∈ D(F ).
• Given µ ∈ D(E) a probability measure on E, M a Markov transition kernel from E to F , f ∈ B(F ) a measurable function on F and

x ∈ E, we adopt the usual following notations:

µ(f ) = f (x)µ(dx), µM (dx ) = E µ(dx)M (x, dx ).
• Given a random variable X on E we denote by LX its law. We have

∀f ∈ B(E), E[f (X)] = LX (f ).
• Given two measures ν and µ on E, ν is absolutely continuous with respect to µ (we write ν µ) if ν(dx) = 0 implies µ(dx) = 0. Moreover if µ is sigma finite, then the Radon-Nikodym derivative of ν with respect to µ exists and is written dν dµ . • For any k ∈ N * , R k is considered as a topological vector space with norm . and the topological structure (R k , O(R k )) is induced by its norm. For each Θ ⊂ R n θ open subset of R n θ , one denotes by VΘ(θ) = {A ∈ O(Θ)|θ ∈ A} the set of all neighborhoods of θ in Θ.

1 Introduction to F-K models and gradient estimation of F-K flows

Feynman-Kac models

In this section we introduce relevant elements of Feynman-Kac theory which will be used all along this paper. For an extensive presentation of theoretical results related to F-K formulae, we refer the reader to the book [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF], and for a good overview of possible applications in the area of particle methods, to the book [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF].

Let (Ep, Ep) p≥0 be a sequence of measurable spaces, and (Xp) p≥0 be a Markov chain from (Ep-1, Ep-1) to (Ep, Ep) with initial measure µ defined on (E0, E0), verifying:

Xp+1 = Fp+1(Xp, Up+1), where Up+1 ∼ ν
where ν ∈ D(E) is a probability measure on the measurable space (U, σ(U )). The transition kernels of the Markov chain (Xp) p≥0 are (Mp) p≥1 , and moreover (Xp) p≥0 is called the canonical Markov chain. Consider a collection of bounded Ep-measurable positive functions Gp ∈ B + (Ep) such that:

∀k ≥ 0, E k p=0 Gp(Xp) > 0.
To a given Feynman-Kac model (µ, (Mp) p≥1 , (Gp) p≥0 ) and a time index n is associated three measures:

• the updated F-K path measure Qn:

Qn(dx0:n) = n p=0 Gp(xp) n p=1 Mp(xp-1, dxp)µ(dx0) E n p=0 Gp(Xp) ,
• the marginal updated F-K measure ηn is the marginalization of Qn:

ηn(dxn) = Qn(E0 × • • • × En-1 × dxn),
• the unnormalized marginal updated F-K measure γn:

γn(dxn) = E 0 ו••×E n-1 n p=0 Gp(xp) n p=1
Mp(xp-1, dxp)µ(dx0).

We immediately make a few elementary remarks. First of all, we have the relation ηn = γn γn [START_REF] Bolic | Resampling algorithms for particle filters: A computational complexity perspective[END_REF] , where 1 denotes here the constant function equals to 1 everywhere. Secondly, as announced in the introduction, the time evolution equation associated to the marginal predicted Feynman-Kac measure (defined by ηn+1 = ηn Mn+1) can be syntactically expressed using the Boltzman-Gibbs measure: ηn+1 = ΨG n+1 ,ηn Mn+1. Finally, given a test function f ∈ B(En), the updated F-K flow ηn(f ) trough f admits a simple expression as a ratio of two expectations along the canonical Markov chain (Xp) p≥0 :

ηn(f ) = E f (Xn) n p=0 Gp(Xp) E n p=0 Gp(Xp) . (6)
We adopt the following convention X0 = F0(X-1, U0), where U0 ∼ ν, and the initial transition kernel verifies M0(x-1, .) = µ.

Interacting particle methods (IPM)

One of the main issue in F-K formulae theory is to compute an empirical measure that approximates the marginal F-K measure ηn by a weighted sum of Dirac measures i wi,nδ ξ i,n . The most famous methods for computing such an approximation are interacting particle methods (IPM) (also called particle filters or sequential monte carlo methods). IPM can be interpreted from many different points of view. First of all, it can be seen as a McKean linearization of the non-linear updating Equation (5) (see [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF] for details). Secondly, one can see IPM as a way to sample from a distribution M g M (g) when knowing only how to sample from and from the Markov transition kernel M (see [START_REF] Douc | Limit theorems for weighted samples with applications to sequential monte carlo methods[END_REF] for details). Thirdly, it can be thought as a voting method that allows to track some signal in the spirit of general bootstrap methods. For our purpose we will think of IPM as stochastic numerical methods that allow to estimate efficiently integrals of the form [START_REF] Crisan | A survey of convergence results on particle filtering for practitioners[END_REF].

A weighted sample {Ξ N p = (w N i,p , ξ N i,p ) 1≤i≤N } N ≥0 is a triangular array of random variables such that w N i,p > 0, i w N i,p = 1 and i (w N i,p ) 2 P → 0 when N → ∞.
IPM propagate a weighted sample along a Feynman-Kac model in three steps:

Prediction. Define ΞN p+1 = ( wN i,p+1 , ξN i,p+1 ) 1≤i≤N with: wN i,p+1 = w N i,p and ξN i,p+1 id ∼ Mp+1(ξ N i,p , •), Updating. Define ΞN p+1 = ( ŵN i,p+1 , ξN i,p+1 ) 1≤i≤N with: ŵN i,p+1 = wN i,p+1 G p+1 ( ξN i,p+1 ) j wN j,p+1 G p+1 ( ξN j,p+1 ) and ξN i,p+1 = ξN i,p+1 , Selection. Ξ N p+1 = S( ΞN p+1 )
where S is a selection operator. The selection operator is designed to reduce the variance of the weights. In literature, selection operators S are often designed such that the weights have no variance after selection and such that the selection step adds no bias:

S(( ŵj, ξj)1≤j≤N ) = 1 N , ξk 1:N E 1 N i f (ξi) | Ξ P = i ŵif ( ξi), (7)
where ki ∈ {1, • • • , N } denotes the i-th selection index used by the selection operator. Many procedures have been built in order to add a minimum of variance in the selection process and to be computationally efficient. The most popular are the multinomial selection, the stratified selection and the residual selection. We refer the reader to [START_REF] Cappé | Comparaison of resampling schemes for particle filtering[END_REF] for a theoretical comparative study of the main selection schemes, and to [START_REF] Bolic | Resampling algorithms for particle filters: A computational complexity perspective[END_REF] for an experimental analysis of them. In the remainder of this paper, we consider a selection scheme that satisfies Equation [START_REF] Douc | Moderate deviation in particle filtering[END_REF].

IPM generate a sequence of these three steps, which may be sketched in the following diagram:

Ξ N 0 -→ . . . -→ Ξ N p Prediction -→ ΞN p+1 Update -→ ΞN p+1 Selection -→ Ξ N p+1 -→ . . . ( 8
)
where ∀i ∈ {1, . . . , N },

ξ N i,0 iid ∼ µ and w N i,0 = 1/N .
Many asymptotic convergence results (when N → ∞) of IPM-based estimate towards ηn can be derived. One can cite for instance: a law a large number, a L q convergence, a central limit theorem, a large and a moderate deviation result. A review of these results is provided in [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF]. The interested reader can find elementary proofs in [START_REF] Crisan | A survey of convergence results on particle filtering for practitioners[END_REF], and some original approaches in [START_REF] Del Moral | Branching and interacting particle systems. approximations of feynman-kac formulae with applications to non-linear filtering[END_REF][START_REF] Douc | Limit theorems for weighted samples with applications to sequential monte carlo methods[END_REF][START_REF] Crisan | A survey of convergence results on particle filtering for practitioners[END_REF][START_REF] Douc | Moderate deviation in particle filtering[END_REF][START_REF] Chopin | Central limit theorem for sequential monte carlo methods and its application to bayesian inference[END_REF]. For our purpose we will only use the following Law of Large Numbers derived from [START_REF] Douc | Limit theorems for weighted samples with applications to sequential monte carlo methods[END_REF]: Theorem 1 (Law of Large Number for Interacting Particle Methods). Using Procedure [START_REF] Douc | Asymptotics of the maximum likelihood estimator for general hidden markov models[END_REF], for any n ≥ 0 we have:

∀f ∈ Bn, 1 N i f (ξ N i,n ) P -→ ηn(f ),
with (Bp) p≥0 defined by induction as follows: B0 = L 1 (η0), and

Bp+1 = {fp+1 ∈ L 1 (ηp+1)|Mp+1(•, Gp+1fp+1) ∈ Bp}.
In the prediction step, one may sample using a sequence of instrumental kernels (Qp) p≥0 satisfying ∀p ≥ 0, ∀xp ∈ Ep, Qp+1(xp, .)

Mp+1(xp, .), instead of sampling using the transition kernels (Mp) p≥0 of the Feynman-Kac flow. If it is the case the weighting step is slightly different because the weights must be multiplied by the Radon-Nikodym derivative dMp dQp . Such importance sampling approaches can lead to a tighter estimate of ( 6), and can be applied to the methods that we will develop in next chapters of this paper. Nonetheless, in order to keep the clarity of the redaction we do not consider importance sampling via an instrumental kernel in the core of the text. We will give a deeper explanation of this remark in Section 3.1.1.

Gradient estimation of F-K flows

Let us now work with a Feynman-Kac model parameterized by a n θdimensional parameter θ ∈ Θ ⊂ R n θ . For simplicity, we omit the θ index, writing Mp for M θ,p , µ for µ θ , and Gp for G θ,p . We are interested in evaluating ∇ηn(f ), the gradient of the updated Feynman-Kac flow with respect to θ for f in a suitable class of functions. Recalling that ηn(f ) = γn(f ) γn(1) , we have:

∇[ηn(f )] = ∇ γn(f ) γn(1) = ∇[γn(f )] γn(1) -ηn(f ) ∇[γn(1)] γn(1) . (9)
Interacting particle methods are well suited to approximate ηn(f ). Consequently the problem is now to estimate ∇[γn(f )] γn(1) .

As we have already noticed in Equation ( 6), γn(f ) can explicitly be expressed in an integral form: γn(f ) = E f (Xn) n p=0 Gp(Xp) . The term ∇γn(f ) γn(1) that remains to be evaluated in Equation ( 9) can be expressed as the gradient of an expectation along the canonical Markov chain:

∇γn(f ) γn(1) = ∇E f (Xn) n p=0 Gp(Xp) γn(1) . ( 10
)
This simple point of view allows us to apply, under quite general assumptions on (Xp) p≥0 , usual methods for gradient of expectations, such as the infinitesimal perturbation method or the score method.

2 Methods for gradient estimation of F-K flows

Let us state the following assumption:

Assumption 1. For any 0 ≤ p ≤ n:

1. the transition function is smooth which means that for any u ∈ U the function (θ, xp-1) → F θ,p (xp-1, u) is differentiable with respect to (θ, xp-1).

the potential function

Gp and the test function f are differentiable with respect to (θ, xp),

3. for any θ ∈ Θ, the random variable f (Xn) verifies:

∃V ∈ VΘ(0), ∃C θ ∈ L 1 (P), ∀h ∈ V, f θ+h (X θ+h,n ) -f θ (X θ,n ) a.s. ≤ h C θ .
We notice that under Assumption 1, for any 0 ≤ p ≤ n, Xp is differentiable with respect to θ. We write ∇Xp its gradient. We also remark that the random variable f (Xp) is in K 1 (Θ, P) (see Appendix for the definition of K 1 (Θ, P)). Moreover, as the functions Gp are all bounded, f (Xn) n p=0 Gp(Xp) remains in K 1 (Θ, P), which is a result of Proposition 2 of Appendix. Therefore we can use Proposition 1 of Appendix to insert the gradient inside the expectation operator and deduce:

∇γn(f ) = ∇E f (Xn) n p=0 Gp(Xp) = E ∇[f (Xn)] n p=0 Gp(Xp) + f (Xn)∇ n p=0 Gp(Xp) = E (∇f (Xn) + f (Xn)∇Xn) n p=0 Gp(Xp) + E f (Xn) n p=0 G p (Xp)∇Xp + ∇Gp(Xp) Gp(Xp) n p=0 Gp(Xp) (11)
We define an augmented canonical Markov chain (Xp, Zp, Rp) p≥0 by the recursive relations:

Xp = Fp(Xp-1, Up), where Up ∼ ν (12) Zp = ∇Fp(Xp-1, Up) + F p (Xp-1, Up)Zp-1 (13)
Rp = Rp-1 + G p (Xp)Zp + ∇Gp(Xp) Gp(Xp) , (14)
where we used the usual matrix notation for the derivatives, ie. Zp is an

nE p × n θ matrix. F p is
∇γn(f ) γn(1) = E (∇f (Xn) + f (Xn)Zn + f (Xn)Rn) n p=0 Gp(Xp) γn(1) = ζn(∇f + A(f ) + B(f )) (15)
where A and B are two linear operators such that A(f )(x, z, r) = f (x)z and B(f )(x, z, r) = f (x)r, and ζn is the marginal Feynman-Kac measure associated to the potential functions (Gp) p≥0 and to the augmented canonical Markov chain (Xp, Zp, Rp) p≥0 .

Plugging together the Equation ( 15) and ( 9) leads to the following estimation of the gradient of Feynman-Kac flow:

∇[ηn(f )] = ∇ γn(f ) γn(1) = ∇[γn(f )] γn(1) -ηn(f ) ∇[γn(1)] γn(1) = ζn(∇f + A(f ) + B(f )) -ηn(f ) ζn(A(1) + B(1)) = ζn(∇f + A(f ) + B(f ) -ηn(f )B(1)). (16)
The right hand side of Equation ( 16) can be evaluated using a numerical method that approximates the marginal Feynman-Kac measure ζn. For example, in the following theorem we use Theorem 1 to estimate it using an Interacting Particle Method: Theorem 2 (Law of Large Number for gradient estimation using IPA).

Under Assumption 1 and for any f ∈ Bn, we have the following asymptotic estimator for the gradient of the Feynman-Kac flow ηn(f ):

1 N N i=1 ∇f (x N i,n ) + f (x N i,n )z N i,n + f (x N i,n ) r N i,n - 1 N j r N j,n P -→ ∇ηn(f ),
where

Ξ N n = ((x N i,n , z N i,n , r N i,n ), w N i,n
) is the weighted sample obtained using an IPM on the Feynman-Kac model associated to the potential functions (Gp) p≥0 and to the augmented canonical Markov chain (Xp, Zp, Rp) p≥0 .

Proof. Starting from Equation ( 16) and using two times Theorem 1 gives:

∇[ηn(f )] = ζn(∇f + A(f ) + B(f ) -ηn(f )B(1)) P ←- 1 N N i=1 ∇f (x N i,n ) + 1 N N i=1 f (x N i,n )z N i,n + 1 N N i=1 f (x N i,n )r N i,n - 1 N N i=1 f (x N i,n ) 1 N j r N j,n = 1 N N i=1 ∇f (x N i,n ) + f (x N i,n )z N i,n + f (x N i,n ) r N i,n - 1 N j r N j,n
The estimator given in Theorem 2 admits a recursive implementation whose computational complexity is linear in the number of particles N . A pseudo-code for its computation is defined at any stage p by: 1. Sample independently ui,p ∼ ν and compute:

xi,p = Fp(xi,p-1, ui,p) zi,p = ∇Fp(xi,p, ui,p) + F p (xi,p, ui,p)zi,p-1 ri,p = ri,p-1 + G p (xi,p)zi,p + ∇Gp(xi,p) Gp(xi,p)
2. Compute the weights ŵi,p = Gp(x i,p ) j Gp(x j,p ) , 3. Define the selection indexes (ki) 1≤i≤N using a consistent selection scheme based on the weights ( ŵi,p)1≤i≤N , and for any i ∈ {1, . . . , N } set (xi,p, zi,p, ri,p) = (x k i ,p , zk i ,p , rk i ,p ).

To study more precisely the asymptotic behavior of the estimator given in Theorem 2, one can apply any of the asymptotic results about Interacting Particle Method (see for example [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF]). For our purpose, we note that this estimator is unbiased and that the asymptotic variance is linked with the variance of the augmented canonical Markov chain (Xp, Zp, Rp) p≥0 (see [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF] or [START_REF] Douc | Limit theorems for weighted samples with applications to sequential monte carlo methods[END_REF] for details). In particular, to be efficient, this algorithm requires that the variances of Zp and Rp do not increase too much as p increases. This is the only possible limitation to the practical efficiency of this algorithm.

The estimator given in Theorem 2, is new in the field of Feynman-Kac models. In [START_REF] Poyiadjis | Particle filter as a controlled markov chain for on-line parameter estimation in general state space models[END_REF] the authors use an Infinitesimal Perturbation Analysis (IPA) to compute the gradient of a filtering flow in the context of HMMs. Nonetheless, they do not apply IPA to the Markov chain (Xp, Zp, Rp) p≥0 but to the empirical Markov chain (Xp, x1:N,p, w1:N,p) p≥0 . Their approach fails to provide an efficient estimator of the Feynman-Kac flow gradient because the variance of their estimator increases as the number of particles N increases (which is not the case for the estimator given in Theorem 2). In fact, in papers [START_REF] Fichoud | Particle-based methods for parameter estimation and tracking : numerical experiments[END_REF][START_REF] Guyader | A particle implementation of the recursive mle for partially observed diffusions[END_REF][START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF] the authors develop algorithms which are very close to the recursive implementation of the estimator given via IPA. Let us now explain precisely the relations.

To compute ∇γn(f ) in Equation ( 11) one could have alternatively used a score method. Under suitable assumptions, one can apply Proposition 3 of Appendix and have:

∇γn(f ) = E ∇f (Xn) + f (Xn) n p=0 ∇mp(Xp-1, Xp) mp(Xp-1, Xp) + ∇Gp(Xp) Gp(Xp) n p=0 Gp(Xp) . ( 17
)
This leads to define a new augmented canonical Markov chain (Xp, Sp) p≥0 by:

Xp = Fp(Xp-1, Up), where Up ∼ ν Sp = Sp-1 + ∇mp(Xp-1, Xp) mp(Xp-1, Xp) + ∇Gp(Xp) Gp(Xp) .
We immediately see that:

∇γn(f ) γn(1) = ςn(∇f + C(f )),
where C is a linear operator verifying C(f )(x, s) = f (x) s, and ςn is the marginal Feynman-Kac measure associated to the potential functions (Gp) p≥0 and to the augmented canonical Markov chain (Xp, Sp) p≥0 . Plugging together the equation above and Equation ( 9) leads to the following estimation of the Feynman-Kac flow gradient:

∇[ηn(f )] = ∇[γn(f )] γn(1) -ηn(f ) ∇[γn(1)] γn(1) = ςn(∇f + C(f ) -ηn(f )C(1)).
Using a particle approximation of ςn immediately leads to the following result:

Theorem 3 (Law of Large Number for gradient estimation using Score).

Under the following assumptions, for any 0 ≤ p ≤ n:

1. the transition kernels Mp(xp-1, dxp) admit a differentiable and bounded density mp(xp-1, xp) with respect to a probability measure ρp ∈ D(Ep),

2. the potential functions Gp are differentiable with respect to θ,

the test function

f : En → R is in K 1 (Θ, ρn) ∩ Bn,
we have:

1 N N i=1 ∇f (x N i,n ) + f (x N i,n ) s N i,n - 1 N N j=1 s N j,n P -→ ∇ηn(f ),
where

Ξ N n = ((x N i,n , s N i,n ), w N i,n
) is the weighted sample obtained using an IPM on the Feynman-Kac model associated to the potential functions (Gp) p≥0 and to the augmented canonical Markov chain (Xp, Sp) p≥0 .

The estimator given in Theorem 3 admits a recursive implementation with computational complexity linear in the number of particles N . A pseudo-code for its computation is defined at any stage p by: 1. Sample independently ui,p ∼ ν and compute:

xi,p = Fp(xi,p-1, ui,p) si,p = si,p-1 + ∇mp(xi,p-1, xi,p) mp(xi,p-1, xi,p) + ∇Gp(xi,p) Gp(xi,p)
2. Compute the weights ŵi,p = Gp(x i,p ) j Gp(x j,p ) , 3. Define the selection indexes (ki) 1≤i≤N using a consistent selection scheme based on the weights ( ŵi,p)1≤i≤N , and for any i ∈ {1, . . . , N } set (xi,p, si,p) = (x k i ,p , sk i ,p ).

This last Algorithm is nearly the same as the one presented in [START_REF] Fichoud | Particle-based methods for parameter estimation and tracking : numerical experiments[END_REF][START_REF] Guyader | A particle implementation of the recursive mle for partially observed diffusions[END_REF][START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF]. The weak point of this approach is that it uses a score method to estimate the sensitivity of the Markov chain (Xp) p≥0 with respect to the parameter θ, and score methods are known to have more variance than IPA methods. A second remark is that our presentation of this algorithm differs from the presentation given in articles [START_REF] Fichoud | Particle-based methods for parameter estimation and tracking : numerical experiments[END_REF][START_REF] Guyader | A particle implementation of the recursive mle for partially observed diffusions[END_REF][START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF]. In particular, the authors of [START_REF] Fichoud | Particle-based methods for parameter estimation and tracking : numerical experiments[END_REF][START_REF] Guyader | A particle implementation of the recursive mle for partially observed diffusions[END_REF][START_REF] Cérou | Stochastic particle methods for linear tangent filtering equations[END_REF] wonder how analyze their algorithm (for example how to obtain a Central Limit Theorem), whereas it is straightforward if one follows the presentation of this article.

3 Numerical experiments and applications

Remarks about practical implementations and possible improvements of the proposed algorithms

In this subsection we discuss some variants and extensions of the Algorithm 2. In a first subsection we explain two variance reduction techniques. Then we explain how to compute differentials of Feynman-Kac flow gradients of orders higher than one. In particular we explicitly show how to adapt the Algorithm 2 to evaluate the Hessian matrix of a Feynman-Kac flow. In a third subsection, we explain how to use the Algorithm 2 if one cannot compute analytically the gradient of functions of the model.

Variance reduction methods

Variance reduction techniques for Interacting Particle Methods have been intensively studied during the last years. See for example the Rao-Blackwellised particle filters [START_REF] Doucet | Raoblackwellised particle filtering for dynamic bayesian networks[END_REF] or the use of auto adaptive minimal variance instrumental kernels [START_REF] Chan | Optimisation of particle filters using simultaneous perturbation stochastic approximation[END_REF]. These methods can be used to accelerate the convergence of the algorithms described in this paper, but their exposition is out the scope of this contribution. In the following, we explain how to use an instrumental kernel and give a slight improvement (in term of variance) of the estimator given in Theorem 2.

To reduce the variance of an estimator obtained using an Interacting Particle Method, one can sample using a family of instrumental transition kernels (Qp) p≥1 instead of sampling from the family of canonical transition kernels (Mp) p≥1 :

x N i,p+1 i.i.d. ∼ Qp+1(x N i,p , .). (18)
Indeed, if the instrumental transition kernels are absolutely continuous with respect to the canonical transition kernels: ∀p ≥ 0, ∀xp ∈ Ep, Qp+1(xp, .) Mp+1(xp, .), then the Radon-Nikodym derivative dMp dQp of (Mp) p≥0 with respect to (Qp) p≥0 is well defined. Under this assumption, the estimator of Theorem 2 and the procedure explained in Section 2 are still valid if one uses the prediction rule given in Equation ( 18 A second way of reducing the variance could be to use a particle approximation available before the selection step. Indeed, the estimator of ∇ηn(f ) given in Theorem 2 is obtained by evaluating the right hand side of Equation ( 16) after the selection step. Nonetheless one can also estimate expectancies by a weighed particle approximation, which means before the selection step. As the selection step is unbiased, using selected particles adds a variance term to the estimator. Consequently, a more accurate estimator of ηn(f ) could be provided by its weighed particle approximation. Therefore, Theorem 1 can be rewritten as

∀f ∈ Bn, i ŵN i,n f ( ξN i,n ) P -→ ηn(f ),
with the same (Bp) p≥0 spaces, and Theorem 2 can also be written with a weighed particle approximation:

Theorem 4 (Law of Large Number for gradient estimation using IPA).

Under Assumption 1 and for any f ∈ Bn, we have the following asymptotic estimator for the gradient of the Feynman-Kac flow ηn(f ):

N i=1 ŵN i,n ∇f (x N i,n ) + f (x N i,n )ẑ N i,n + f (x N i,n ) rN i,n - j ŵN j,n rN j,n P -→ ∇ηn(f ),
where ΞN n = ((x N i,n , ẑN i,n , rN i,n ), ŵN i,n ) is the weighted sample obtained using an IPM on the Feynman-Kac model associated to the potential functions (Gp) p≥0 and to the augmented canonical Markov chain (Xp, Zp, Rp) p≥0 .

Proof. Starting from Equation ( 16) and using two times the weighted version of Theorem 1 (see above) gives:

∇[ηn(f )] = ζn(∇f + A(f ) + B(f ) -ηn(f )B(1)) P ←- N i=1 ŵN i,n ∇f (x N i,n ) + N i=1 ŵN i,n f (x N i,n )ẑ N i,n + N i=1 ŵN i,n f (x N i,n )r N i,n - N i=1 ŵN i,n f (x N i,n ) j ŵN j,n rN j,n = N i=1 ŵN i,n ∇f (x N i,n ) + f (x N i,n )ẑ N i,n + f (x N i,n ) rN i,n - j ŵN j,n rN j,n
The IPA Algorithm is not deeply modified. One has to use predicted particles and weights instead of using selected particles and all the same weights equal to 1 N .

3.1.2 Estimation of higher order differentials. A case study with the Hessian matrix.

The method developed in Section 2 allows to estimate ∇ηn(f ). One should be interested in estimating higher order differentials such as the Hessian matrix H[ηn(f )] of the Feynman-Kac flow. This is an interesting issue in optimization context, because it allows the use of second order optimization methods, as the Newton-Raphson algorithm, that achieve a quadratic rate of convergence to a locally optimal solution. The Algorithm 2 can be adapted to compute differential of ηn(f ) of any order. As an illustration, we explain in the following how to modify the Algorithm 2 in order to approximate the Hessian matrix H[ηn(f )].

The Hessian matrix H[ηn(f )] is a n θ × n θ matrix whose element at the position i, j is ∂ 2 ηn(f )

∂ θ i ∂ θ j
. To simplify the explanation, we suppose in the following that the test function f doesn't depend on the parameter θ, n θ = 1 and for any p ≥ 0, nE p = 1. We also suppose that all the quantities of interest are sufficiently smooth. Under these assumptions, we have:

H[ηn(f )] = ∇ [∇ηn(f )] = ∇ ∇[γn(f )] γn(1) -ηn(f ) ∇[γn(1)] γn(1) = H[γn(f )] γn(1) - ∇γn(f ) γn(1)
∇γn( 1) γn( 1) -∇ηn(f ) ∇γn( 1) γn( 1)

-ηn(f ) H[γn(1)] γn(1) - ∇γn(1) γn(1) 2 = H[γn(f )] γn(1) -2 ∇γn(f ) γn(1)
∇γn( 1) γn( 1)

-ηn(f ) H[γn(1)] γn(1) -2 ∇γn(1)
γn( 1)

2 ( 19
)
To evaluate the right hand side of Equation ( 19), we have to evaluate H[γn(f )]. For this purpose, we can apply two times Proposition 1:

H[γn(f )] = ∇ ∇E f (Xn) n p=0 Gp(Xp) = ∇E (f (Xn)∇Xn + f (Xn)Rn) n p=0 Gp(Xp) = E ∇(f (Xn)∇Xn + f (Xn)Rn) + (f (Xn)∇Xn + f (Xn)Rn)Rn n p=0 Gp(Xp) = E f (Xn)∇X 2 n + f (Xn)(H[Xn] + 2∇XnRn) + f (Xn)(∇Rn + R 2 n ) n p=0
Gp(Xp)

The above relation can be seen as an expectation along the Markov chain (Xp, Zp, Hp, Rp, Tp) p≥0 where the recurrence relations verified by Xp, Zp and Rp have already been explained in Section 2 and the recurrence relations verified by Hp is given by: [START_REF] Poyiadjis | Particle Method for Parameter Estimation in General State Space Models[END_REF] and by Tp verifies:

Hp = H[Xp] = ∇[∇Xp] = ∇[∇Fp(Xp-1, Up) + F p (Xp-1, Up)Zp-1] = H[Fp](Xp-1, Up) + 2∇F p (Xp-1, Up)Zp-1 +F p (Xp-1, Up)Hp-1 + F p (Xp-1, Up)Z 2 p-1
Tp = ∇Rp = ∇ Rp-1 + G p (Xp) Gp(Xp) Zp + ∇Gp(Xp) Gp(Xp) = Tp-1 + ∇ G p (Xp) Gp(Xp) Zp + G p (Xp) Gp(Xp) Hp] + ∇ ∇Gp(Xp) Gp(Xp) = Tp-1 + ∇[G p (Xp)] Gp(Xp) Zp - G p (Xp)∇[Gp(Xp)] Gp(Xp) 2 Zp + G p (Xp) Gp(Xp) Hp + ∇[∇Gp(Xp)] Gp(Xp) - ∇Gn(Xp)∇[Gp(Xp)] Gp(Xp) 2 = Tp-1 + ∇G p (Xp) + G p (Xn)Zp Gp(Xp) Zp - G p (Xp)∇Gp(Xp) + G p (Xp) 2 Zp Gp(Xp) 2 Zp + G p (Xp) Gp(Xp) Hp + H[Gp](Xp) + (∇Gp) (Xp)Zp Gp(Xp) - ∇Gp(Xp) 2 + ∇Gp(Xp)G p (Xp)Zp Gp(Xp) 2 (21)
Using an Interacting Particle Method on the Feynman-Kac model associated to the Markov chain (Xp, Zp, Hp, Rp, Tp) p≥0 and to the potential functions (Gp) p≥0 allows to numerically estimate the right hand side of Equation [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF]. A sketch of the pseudo code is:

1. Sample independently ui,p ∼ ν and compute ξ1:N,p = (x1:N,p, z1:N,p, h1:N,p, r1:N,p, t1:N,p) using the recurrence relations ( 12),( 13),( 14),( 20) and ( 21), 2. Compute the weights ŵ1:N,p = Gp(x 1:N,p ) j Gp(x j,p ) , 3. Define the selection indexes (ki) 1≤i≤N using a consistent selection scheme based on the weights ( ŵi,p)1≤i≤N , and for any i ∈ {1, . . . , N } sets ξ1:N,p = ξk 1:N ,p .

What to do if one cannot compute analytically gradients of functions needed for algorithms?

differentials, it is possible to use finite difference estimation. For example this is the case when transition functions are the result of a full computer program. This situation is quite frequent in industrial applications.

In such a situation the Algorithm 2 still works, but we need to evaluate ∇Fp or F p by a finite difference method:

∇Fp ≈ F θ+h i ,p -F θ-h i ,p 2h 1≤i≤n Θ
where hi is a vector for which all coordinates are null except the i-th which is equal to 0 < h 1. Note that this evaluation requires the computation of 2nΘ transition evaluations and this could be too computationally expensive in situations where nΘ 1 and where sampling a transition is time consuming.

Application to gradient evaluation and parameter identification in HMM

In this section we briefly recall the notations and definitions used in Section 1.1 of this paper. We show how to apply the estimators given in Theorem 2 to the problem of likelihood evaluation and maximization in HMMs.

A HMM is built from a state process and an observation process. The state process (Xp) p≥0 is an homogeneous Markov chain with initial probability measure µ(dx0) and Markov transition kernel K(dxp+1|xp). The observation process (Yp) p≥1 is linked with the state process by the conditional probability measure P(Yp ∈ dyp|Xp = xp) = g(yp, xp)λ(dyp) and is conditionally independent given the state process: i = j ⇒ P(Yi ∈ dyi, Yj ∈ dyj|Xi = xi, Xj = xj) = P(Yi ∈ dyi|Xi = xi)P(Yj ∈ dyj|Xj = xj). Moreover, we consider HMMs for which the state and the observation processes are parameterized by θ ∈ Θ, where Θ ⊂ R n θ is an open subset of R n θ . We present two classical examples of parameterized HMMs. They are frequently used to compare performances of algorithms dealing with identification and filtering in HMMs.

Example 3.1 (Autoregressive model). The autoregressive model of order 1 (AR1) is the simplest representation of the class of linear-Gaussian

HMMs:

Xp = φXp-1 + σUp, X0 ∼ N 0, σ 2 1 -φ 2 Yp = ρXp + βVp, (22)
where Up i.i.d.

∼ N (0, 1) and Vp i.i.d.

∼ N (0, 1) are two mutually independent and independent of the initial state X0 sequences of random variables. In this example, θ = (φ, σ, ρ, β) is a four-dimensional parameter and Θ = (R+) 4 . For an AR model the filtering density (and therefore the likelihood function and its gradient) can be computed exactly using a Kalman filter. This model will allow us to measure the bias of algorithms given in Section 2, by comparing them with the true estimation.

Example 3.2 (Stochastic Volatility). Hull and al. generalized the Black-Scholes option pricing formula to allow for stochastic volatility. Their formula has emerged as the dominant approach. In fact it is a typical nonlinear and nongaussian HMM. The dynamics are:

Xp = φXp-1 + σUp, X0 ∼ N 0, σ 2 1 -φ 2 Yp = β exp Xp 2 Vp, (23)
where Up

i.i.d.
∼ N (0, 1) and Vp

i.i.d.
∼ N (0, 1) are two mutually independent and independent of the initial state X0 sequences of random variables. It deals with a three-dimensional parameter θ = (φ, σ, β) and Θ = R 3 + . This model will allow us to assess the capacity of the Algorithm 2 to deal with non linear and non gaussian dynamics.

Let us now suppose we are given a sequence of successive realizations of the observation process (y p,θ * ) p≥1 that were obtained using an unknown parameter denoted θ * . In the following, we make a slight abuse of notation by denoting g p,θ (xp) = g θ (y p,θ * , xp). We recall that the objective of parameter identification in HMMs is to recover θ * using the sequence of observations (y p,θ * ) p≥1 . We have seen in the introduction that the parameter estimation problem can be turned into maximizing the log-likelihood ln. In practice, this can be achieved with a stochastic gradient ascent with respect to θ on ln. We also recall that the gradient of the log-likelihood function may be written as:

∇ln(θ) = n p=1 ∇[π p|p-1,θ (g p,θ )] π p|p-1,θ (g p,θ ) . ( 24
)
We have also seen in Section 1.1 that a HMM can be viewed as a particular case of Feynman-Kac model. Indeed, taking for any p ≥ 0, Mp = K and Gp = gp, implies that ηp+1 = π p+1|p is the marginal predicted measure of the Feynman-Kac model associated to (Gp, Mp) p≥0 . Consequently we have:

∇[π p|p-1 (g p,θ )] = ∇ηp(Gp). (25)
Theorem 2 has been designed to evaluate the gradient of the updated Feynman-Kac flow ∇ ηp(f ) and not the gradient of the predicted Feynman-Kac flow ∇ηp(f ). Nonetheless, Theorem 2 can be easily adapted to evaluate the gradient of the predicted Feynman-Kac flow ∇ηp(f ). Indeed the algorithm remains the same, but the evaluation is done using the predicted particles system:

1 N N i=1 ∇f (x N i,p ) + f (x N i,p )z N i,p + f (x N i,p ) r N i,p-1 - 1 N j r N j,p-1 P -→ ∇ηn(f ). ( 26
)
To evaluate the right hand side of Equation (25), we can apply the evaluation provided in Equation ( 26) taking f = Gp, this leads to the following estimator:

1

N N i=1 ∇Gp(x N i,p ) + G p (x N i,p )z N i,p + Gp(x N i,p ) r N i,p-1 - 1 N j r N j,p-1
.

Consequently, we have:

N i=1 ∇Gp(x N i,p ) + G p (x N i,p )z N i,p + Gp(x N i,n ) rN i,p -1 N j rN j,p j Gp(x N j,p ) P -→ ∇[π p|p-1,θ (g p,θ )] π p|p-1,θ (g p,θ ) . (27)
The above procedure shows how to estimate the gradient of the loglikelihood ∇ln. If we want to recover the unknown parameter θ * , we have seen that it is possible to compute a stochastic gradient ascent using ∇ln. A pseudo-code of the application of Algorithm 2 to the parameter estimation in HMMs when n observations are available, is given through the repetition for increasing time step k of these two points:

1. estimate ∇ln(θ k ) using Equation 27 and Equation 242. update the parameter

θ k+1 = θ k + k ∇ln(θ k )
where k k = ∞ and k 2 k < ∞. A main advantage of this algorithm is that it can be used on-line.

Numerical validation on an AR 1 model

In this section we present results of numerical experiments on a gradient estimation problem of the log-likelihood associated to the linear gaussian Hidden Markov Model defined in Example 3.1 with θ * = (φ * , σ * , ρ * , β * ) = (0.8, 0.5, 1.0, 1.0). The Kalman filter allows to compute the exact value of the filtering distribution, and therefore the exact value of the log-likelihood function ln(θ) and its gradient ∇ln(θ) for any values of the parameter of the model θ. We know that for linear gaussian models the Algorithms obtained from the two Theorems 2 and 3 have no interest because the gradient estimation problem is exactly solved by using the Kalman filter. Nonetheless, as the exact value of the log-likelihood gradient ∇ln(θ) is known, it allows us to evaluate the mean quadratic error e N n of the approximation of ∇ln(θ) by the value ∇ln N (θ) computed with Algorithms 2 and 3 using N particles.

During numerical experiments, we evaluate separately the bias and the variance of the gradient estimator ∇ln In the following three subsections, we evaluate the performances of Algorithms 2,3 and [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF]. The experiments were realized with a 1.6 Ghz PC under the software MatLab. We have not optimized the implementation of the three algorithms, in particular we have not implemented the acceleration method proposed in [START_REF] Doucet | Toward practicle n 2 monte carlo: the marginal particle filter[END_REF] for the marginal particle filter as used in [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF]. The acceleration method adds a numerical approximation through the use of a tree-based numerical method to compute quickly the linear combination of the kernel evaluation. Moreover we study the behavior of Algorithms 2 and 3 along the two variables N and n: when the number of particle N tends to infinity and when the number of times step n tends to infinity. [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] In this section we present a comparative study of performances of Algorithms 2, 3 and [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] in terms of bias, variance and CPU time. We evaluate the log-likelihood at time step n = 50 and at point θ = (φ, σ, ρ, β) = (0.7, 0.4, 0.9, 0.9). For each algorithm, we compute m = 500 log-likelihood estimates ( ∇l 50,k N (θ)) 1≤k≤500 , and evaluate empirically the bias and the variance using Equation (28). For Algorithms 2 and 3, we use successively N = 5 10 2 particles and N = 1 10 4 particles, and for the Algorithm [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] we only use N = 5 10 2 particles as it is practically impossible to apply it using more than N = 10 3 particles. We also measure the computational time (in second) needed to evaluate one log-likelihood ∇l 50,k N (θ). The results are reported in Table 1.

Comparisons with the state of art method

As it is expected, the first remark is that the bias is null. In these numerical simulations it is not exactly null because we evaluate it using a finite number of simulations m = 500. But we see that the bias is always inferior to the standard deviation in one or two orders of dimension. A second conclusion is that the all three Algorithms 2, 3 and [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] achieve a rather small quadratic error e 0.5 < 10 -2 , thus they are likely to be asymptotically consistent. A third conclusion is that the method [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] 8.8 10 -3 6.0 10 -2 X CPU time (s.) 5 10 2 1.8 10 -1 1.9 10 -1 21 10 4

2.1 2.3 X bias 5 10 2 1.3 10 -3 4.0 10 -3 1.0 10 -4 10 4 2.0 10 -4 1.1 10 -3 X standard deviation 5 10 2 2.3 10 -2 6.6 10 -2 2.2 10 -2 σ 10 4 7.9 10 -3 2.0 10 -2 X CPU time (s.) 5 10 2 1.7 10 -1 1.9 10 -1 23 10 4 2.1 2.3 X bias 5 10 2 6.0 10 -4 1.2 10 -3 8.0 10 -4 10 4 3.0 10 -4 1.0 10 -4 X ρ standard deviation 5 10 2 1.5 10 -2 1.5 10 -2 6.9 10 -2 10 4 6.0 10 -3 5.7 10 -3 X CPU time (s.) 5 10 2 1.9 10 -2 2. Indeed, even for a small number of particles N = 500, the CPU time needed to compute a log-likelihood approximation is 10 2 -times the CPU time needed for Algorithms 2, 3. The fourth conclusion is that for a fixed number of particles, the method [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF] performs slightly better than the two others algorithms. However, as we can notice in experiments using N = 10 4 particles, for a fixed CPU time the two others algorithms strongly outperform the Algorithm [START_REF] Poyadjis | Particle methods for optimal filter derivative: Application to parameter estimation[END_REF].

Estimation behavior when the number of particles increases

We have seen in Subsection 3.3.1 that the estimators of Algorithms 2 and 3 achieve a rather small quadratic error for N = 10 4 particles. Moreover, we have noticed in Section 2 that it is possible to deduce a Central Limit Theorem for these estimators. It is interesting to study experimentally the dependance of the quadratic error e N n in the number of particles N . To evaluate the quadratic error of each algorithm, we The main conclusion is that the experimental results strengthen the consistence result given in Section 2: the error tends to zero as the number of particle tends to infinity. A second remark is that the two methods have the same performances for coefficients in the observable equation. In fact, this is obvious because the algorithms are the same for coefficients belonging only to the potentials functions (Gp) p≥0 . A third remark is that the IPA estimator (Theorem 2) performs slightly better than the Score estimator (Theorem 3) for coefficients in the state equation (in F-K context, it means coefficients belonging to the transition kernels (Mp) p≥0 ). In this example it is particulary true for the coefficient σ which is the noise of the state equation.

Estimation behavior when the number of time steps increases

IPA

Score n 10 10 2 10 3 10 10 2 10 3 bias 2.0 10 -4 1.7 10 -3 6.0 10 -4 3.0 10 -4 4.0 10 -4 2.3 10 -3 φ standard deviation 1.8 10 -2 3.3 10 -2 2.7 10 -2 2.7 10 -2 2.6 10 -2 2.5 10 -2 CPU time (s) 5.0 10 -2 5.0 10 -1 5.2 5.0 10 -2 5.5 10 -1 5.3 bias 9.0 10 -4 2.5 10 -3 1.3 10 -3 4.7 10 -3 7.0 10 -4 1.2 10 -3 σ standard deviation 2.3 10 -2 2.2 10 -2 2.2 10 -2 6.5 10 -2 6.1 10 -2 5.6 10 -2 CPU time (s) 4.0 10 -2 5.0 10 -1 5.2 6.0 10 -2 5.3 10 -1 5.6 bias 7.0 10 -4 3.0 10 -4 4.0 10 -4 8.0 10 -4 6.0 10 -4 7.0 10 -4 ρ standard deviation 1.5 10 -2 9.7 10 -3 1.2 10 -2 1.6 10 -2 8.9 10 -3 1.3 10 -2 CPU time (s) 5.0 10 -2 5.3 10 -1 5.2 5.0 10 -2 5.8 10 -1 5.5 bias 7.0 10 -4 1.9 10 -3 2.6 10 -3 8.4 10 -6 2.3 10 -3 8.0 10 -4 β standard deviation 1.7 10 -2 1.9 10 -2 2.3 10 -2 1.8 10 -2 1.9 10 -2 2.0 10 -2 CPU time (s) 6.0 10 -2 5.3 10 -1 5.7 5.0 10 -2 5.8 10 -1 5.9

Table 2: Asymptotic behavior along the number of time steps n of Algorithms 2 and 3 by computing the log-likelihood derivative of Example 3.1. We use N = 10 3 particles. The values of the true parameters are (0.8, 0.5, 1.0, 1.0), and the log-likelihood derivatives are computed at (0.7, 0.4, 0.9, 0.9).

One of the main issue in Feynman-Kac flow estimation is the stability of the numerical approximation when the number of time steps n increases. In the HMM framework this corresponds to an increasing number of observations y1:n. To assess the stability of the estimator of Algorithms 2 and 3, we evaluate the quadratic error for each algorithm by computing m = 500 log-likelihood estimates ( ∇l 50,k N (θ)) 1≤k≤500 at point θ = (φ, σ, ρ, β) = (0.7, 0.4, 0.9, 0.9) with N = 10 3 particles and for different values of the number of time steps n ∈ {10, 10 2 , 10 3 }. The results are presented in Table 2.

The conclusion of this last study is that the quadratic error e N n is independent of the number of time steps n, which experimentally confirms the stability of the two Algorithms 2 and 3. In this subsection, we are interested in evaluating the log-likelihood function for the stochastic volatility model presented in Example 3.2. As this model is neither linear nor gaussian, the log-likelihood function cannot be evaluated using a Kalman filter. We have seen that parameter identification can be reduced to searching the zero of the log-likelihood gradient function. Consequently, we have plotted the value of one evaluation of ∇ln N (θ) for n = 5 10 3 , for different values of parameter θ, with θ * = (φ , σ , β ) = (0.8, 0.5, 1), and using N = 5 10 4 particles. Results are represented in Figure 2.

Numerical experiments on a stochastic volatility model

These plots show that the estimation is a little bit noisy but the estimated derivatives are equal to zero at a point belonging to a neighborhood of the true parameter, which confirms that the derivative of log-likelihood function is well evaluated.

Appendix

Let (Ω, F , P) be a probability space. For any k ∈ N * , R k is considered as a topological vector space with norm and the topological structure (R k , O(R k )) induced by its norm. In the following {(ω, θ) → X θ (ω)} is a family of R m -valued random variables indexed by a parameter θ ∈ Θ, where Θ ∈ O(R n θ ) is an open subset of R n θ . For any θ ∈ Θ, one denotes by VΘ(θ) = {A ∈ O(Θ)|θ ∈ A} the set of all neighborhoods of θ in Θ.

We adopt the following notation: Notation 1. K q (Θ, P) denotes the set of all functions {(ω, θ) → X θ (ω)} such that for any θ ∈ Θ:

1. the function {ω → X θ (ω)} is measurable, 2. for almost all ω ∈ Ω the function {θ → X θ (ω)} is differentiable at θ.

We write ∇X θ its gradient, 3. ∃V ∈ VΘ(0), ∃C θ ∈ L q (P), ∀h ∈ V, X θ+h -X θ a.s.

≤ h C θ .

Proposition 1 (Dominated convergence for gradient).

If X ∈ K 1 (Θ, P) then for any θ ∈ Θ, the function {ω → ∇ θ X(ω)} is in L 1 (P), the function {θ → E[X θ ]} is differentiable and verifies ∇E[X θ ] = E[∇X θ ].

Proof. Given a sequence hn of vectors, we set F θ,n (ω) =

X θ+hn (ω)-X θ (ω) hn . For any sequence hn → 0, we have:

∇E[X θ ] = lim n E[X θ+hn ] -E[X θ ] hn = lim n E X θ+hn -X θ hn = lim n E[F θ,n ].
For almost all ω ∈ Ω, F θ,n (ω) ≤ g θ (ω), using the dominated convergence theorem we have limn The following assertions are true:

1. For any q ∈ N ∪ {∞}, a differential and locally contracting function is in K q (Θ, P). In particular C 1 (Θ) ⊂ K q (Θ, P).

2. For any q ∈ N ∪ {∞}, K q (Θ, P) is a vectorial space.

3. For any q ∈ N ∪ {∞}, if X ∈ K q (Θ, P) and (Y θ ) θ∈Θ is locally bounded:

∀θ ∈ Θ, ∃B θ > 0, ∃V ∈ VΘ(0), ∀h ∈ V, Y θ+h ∞ < B θ , then XY ∈ K q (Θ, P).

Proof. 1. Given a family of function (X θ ) θ∈Θ differentiable with respect to θ and locally contracting with local contraction constant C θ , for any q > 0 we have C θ (ω) = C θ ∈ L q (P).

As any continuously differentiable function taking its values in a locally compact vectorial space is locally contracting, the conclusion follows from the remark above.

2. For any q > 0, take λ ∈ R, X ∈ K q (Θ, P) and Y ∈ K q (Θ, P), for almost all ω ∈ Ω we have: As L q (P) is a vectorial space (|λ|C X θ (ω) + C Y θ (ω)) ∈ L q (P).

3. For almost all ω ∈ Ω we have:

X θ+h (ω)Y θ+h (ω) -X θ (ω)Y θ (ω) ≤ X θ+h (ω) -X θ (ω) B θ ≤ h B Y θ C X θ (

  ) and the following updating rule: ŵi,p = dMp dQp (xi,p-1, xi,p)Gp(xi,p) j dMp dQp (xj,p-1, xj,p)Gp(xj,p) .

N 2 = 2 ,

 22 (θ) by computing independently m-times the quantity ∇ln N (θ) and by using the following bias/variance decomposition: e N n (θ) = E ∇ln(θ) -∇ln N (θ)) E ∇ln(θ) -∇ln N (28) where E ∇ln(θ) -∇ln N (θ) is called the bias of the estimator and var ∇ln N (θ) the variance.

  is very time consuming when one doesn't use any acceleration technique.

Figure 1 :

 1 Figure 1: Asymptotic behavior along the number of particles N of Algorithms 2 and 3 by computing the log-likelihood derivative of Example 3.1. We use 10 2 , 10 3 , 10 4 , and 10 5 particles, and 50 time steps. The values of the true parameters are (0.8, 0.5, 1.0, 1.0), and the log-likelihood derivatives are computed at (0.7, 0.4, 0.9, 0.9). Left up: φ. Right up: σ. Left down : ρ. Right down : β.

Figure 2 :

 2 Figure 2: Graphical representation, for some values of θ, of the evaluation of the log-likelihood gradient function ∇l n N (θ) associated to the stochastic volatility model presented in Example 3.2 at time step n = 5 10 3 and using the estimator given in Algorithm 2. The true parameter is θ * = (φ , σ , β ) = (0.8, 0.5, 1) and we use N = 5 10 4 particles. Left up: φ. Right up: σ. Down : β.

  E[F θ,n ] = E[limn F θ,n ]. And, by differentiability of the function {θ → X θ (ω)}, we have limn F θ,n a.s.= ∇X θ . Example 4.1. Take Ω = [0, 1], P = U([0, 1]), Θ =]0, 1[ and X θ (ω) = 1 [θ,1] . For any θ ∈ Θ and any ω = θ, the function {ω → X θ (ω)} is differentiable with value ∇X θ (ω) = 0. Note that X θ+h (ω) -X θ (ω) = 1 [θ+h,1] -1 [θ,1] = -1 [θ,θ+h[ implies sup h |X θ+h (ω)-X θ (ω)| |h| =1|ω-θ| , consequently X is not in K 1 (Θ, P).

Example 4 . 2 .

 42 Using the same spaces, now takeX θ (ω) = ω-θ b 1 [θ,θ+b[ (ω)+ 1 [θ+b,1[ (ω). X is in K 1 (Θ, P) because |X θ+h (ω) -X θ (ω)| ≤ |h| b .Proposition 2 (Stability of K q (Θ, P) spaces).

  λX θ+h (ω) + Y θ+h (ω) -λX θ (ω) -Y θ (ω) ≤ |λ| X θ+h (ω) -X θ (ω) + Y θ+h (ω) -Y θ (ω) ≤ (|λ|C X θ (ω) + C Y θ (ω)) h .

ω) Proposition 3 ( 1 ,

 31 Score Method applied to a Markov chain).Let n ∈ N be an integer, (Xp) p≥0 be a Ep-valued Markov chain and f a function f : × n p=0 Ep → R. If for any 0 ≤ p ≤ n, the transition kernels Kp(xp-1, dxp) = k(xp-1, xp)νp(dxp) admit a differentiable and bounded density relative to a probability measure νp ∈ D(Ep), and f ∈ K 1 (Θ, ∧ n p=0 νp), then:∇E[f (Xn, . . . , X0)] = E f (Xn, . . . , X0) Xp) .Proof. First we note that E[f (Xn, . . . , X0)] = f (xn, . . . , x0) n p=0 Kp(xp-1, dxp). As Kp(xp-1, dxp) = k(xp-1, xp)νp(dxp), we have f (xn, . . . , x0) n p=0 Kp(xp-1, dxp) = f (xn, . . . , x0) n p=0 kp(xp-1, xp)νp(dxp). The functions kp(xp-1, xp) are bounded and differentiable, and f ∈ K 1 (Θ, ∧ n p=0 νp), then thanks to the result 3 of Proposition 2, we have: f n p=0 kp ∈ K 1 (Θ, ∧ n p=0 νp). Consequently, using Proposition 1, we have: ∇ f (xn, . . . , x0) n p=0 kp(xp-1, xp)νp(dxp) = f (xn, . . . , x0)∇[ n p=0 kp(xp-1, xp)] n p=0 νp(dxp) = f (xn, . . . , x0) n p=0 ∇kp(x p-1 ,xp) kp(x p-1 ,xp) n p=0 kp(xp-1, xp)νp(dxp) = E f (Xn, . . . , X0) n p=0 ∇kp kp (Xp-1, Xp) Remarking (k(xp-1, xp) = 0) ⇒ (∇k(xp-1, xp) = 0), we see that the ratio ∇kp(x p-1 ,xp) kp(x p-1 ,xp) is well defined.

Table 1 :

 1 Results computed at φ = 0.7 whereas φ = 0.8 with 500 simulations, 500 particles and 50 time steps.

				1 10 -1	19
		10 4	2.2	2.3	X
	bias	5 10 2 1.7 10 -4 3.9 10 -3 4.5 10 -3
		10 4	2.0 10 -4 1.0 10 -4	X
	standard deviation 5 10 2 4.4 10 -2 4.3 10 -2 2.7 10 -2
	ρ	10 4	6.2 10 -3 5.7 10 -3	X
	CPU time (s.)	5 10 2 2.0 10 -1 2.1 10 -1	20
		10 4	2.3	2.5	X

* CMAP, Ecole Polytechnique (France), Email: name@cmapx.polytechnique.fr † SequeL team, INRIA futurs Lille (France), Email : firstname.name@inria.fr

If transition functions are smooth (they verify Assumption 1) but one cannot (or doesn't want to) compute analytically the expression of their