Tomas G Moreira
email: tgmoreira@inf.ufrgs.br

Marco A Wehrmeister
email: marcow@joinville.udesc.br

Carlos Eduardo Pereira
email: cpereira@ece.ufrgs.br

Jean-François Pétin
email: jean-francois.petin@cran.uhp-nancy.fr

Eric Levrat
email: eric.levrat@cran.uhp-nancy.fr

Tomás G Moreira

Generating VHDL source code from UML models of embedded systems

Keywords: Embedded systems, system engineering, intelligent maintenance, UML specification, VHDL code generation

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Embedded systems are dedicated system designed to perform a small number of functions. It contains predominantly digital components, consisting in a hardware platform upon which software application execute [START_REF] Micheli | Hardware/Software Co-design[END_REF]. Embedded systems' functionalities can be distributed over different processing nodes. Distributed Embedded Systems (DES) rely on a communication infrastructure constrained by requirements/constraints of embedded systems domain, e.g. timing and energy consumption requirements.

In the industrial domain, DES may be composed by several intelligent components, which make decisions and perform their activities autonomously [START_REF] Pereira | Distributed real-time embedded systems: Recent advances, future trends and their impact in manufacturing plant control[END_REF]. Industrial DES support conventional or innovative functions. The former concerns to simple control functions, whereas the later represents more elaborated functions, e.g. maintenance and prognostic systems' functions to perform Condition Monitoring (CM), Health Assessment (HA), Prognostics (PR), etc. Components' intelligence level is defined by the amount of different services required by the end-user that are implemented as component functions [START_REF] Iung | Engineering process of integrated-distributed shop floor architecture based on interoperable field components[END_REF], [START_REF] Pétin | Distributed intelligent actuation and measurement (IAM) system within an integrated shop-floor organization[END_REF]. As machines do not suddenly fail, they usually pass through a measurable process of degradation before failing. Intelligent maintenance or prognostics systems use information provided by sensors and computing components embedded into equipments. Then, algorithms for health estimation and failure prediction are applied to assess machines' degradation level. Hence, embedded sensors, intelligent actuators and processing elements play a fundamental role in the development of intelligent maintenance systems.

The complexity of this scenario demands new tools and techniques. Increasing the design's abstraction level by using, for instance Model-Driven Engineering (MDE) [START_REF] Mellor | Guest Editors' Introduction: Model-Driven Development[END_REF] techniques, is an interesting approach to deal with the mentioned issues [START_REF] Sangiovanni-Vincentelli | The tides of EDA[END_REF]. Standard graphical languages, e.g. the Unified Modeling Language (UML) 1 , must be used to facilitate the communication of design's intentions to different teams, i.e. software and hardware teams. UML is a high-level design language and is broad enough in scope to model DES. Usually, UML-based MDE approach focus only in the software part of embedded systems by defining a mapping between high-level specification's construction down to software construction using programming languages (e.g. C/C++, Java, etc.). There are many academic works and commercial tools (e.g. Rational Rose 2 and Artisan Studio 3 ) that can generate software code from UML models. However, considering embedded systems' hardware part, only few works address the use of UML to produce Hardware Description Language (HDL) descriptions, as in [START_REF] Mcumber | UML-Based Analysis of Embedded Systems Using a Mapping to VHDL[END_REF], [START_REF] Coyle | From UML to HDL: a Model-Driven Architectural approach to hardware-software co-design[END_REF]. In this sense, the transformation of UML models into HDL code, e.g. using VHDL (Very High Speed Integrated Circuit HDL), is not yet well diffused, opening room for research on this subject.

This works presents an extension to our aspect-oriented MDE approach for DES named Aspect-Oriented Model-Driven Engineering for Real-Time systems (AMoDE-RT) [START_REF] Wehrmeister | An Aspect-Oriented Model-Driven Engineering Approach for Distributed Embedded Real-Time Systems[END_REF]. This work's main contribution is to support automatic generation of VHDL descriptions from UML models. In other words, this work extends AMoDE-RT's supporting tool GenERTiCA (Generation of Embedded Real-Time Code based on Aspects) [START_REF] Wehrmeister | GenERTiCA: A Tool for Code Generation and Aspects Weaving[END_REF], aiming at automatic generation of VHDL descriptions from UML models. The generated VHDL code is intended to be used in FPGA (Field Programmable Gate Array) systems. Thus, the proposed approach allows the behavior of the required system to be described (modeled) and verified (simulated) before synthesis tools translate the design into real hardware (gates and wires).

Additionally, this work proposes an engineering process, which covers from requirements analysis and UML modeling phases to VHDL code generation. The focus is to generate VHDL source code for the logical functions of an embedded system, which so far has only been implemented in software. To validate this work, this paper presents a use case focused on a distributed embedded system (i.e. DES) used for maintenance systems (intelligent components), which integrates both conventional and innovative functions.

This paper is organized as follows: Section 2 provides an overview of the proposed approach to map UML specifications into VHDL source code, whereas Section 3 presents the developed mapping rules. The case study of a valve component system is presented in Section 4. Section 5 provides a review of works related. Finally, conclusion and directions for future work are presented in Section 6

2

Overview of the Proposed Approach

The proposed approach follows the flow proposed by the Aspect-oriented Model-Driven Engineering for Real-Time systems (AMoDE-RT) [START_REF] Wehrmeister | An Aspect-Oriented Model-Driven Engineering Approach for Distributed Embedded Real-Time Systems[END_REF], [START_REF] Wehrmeister | An Aspect-Oriented Approach for Dealing with Non-Functional Requirements in Model-Driven Development of Distributed Embedded Real-Time Systems[END_REF] (see Fig. 1), which uses MDE techniques combined with AO concepts to design DES. AMoDE-RT is supported by GenERTiCA [START_REF] Wehrmeister | GenERTiCA: A Tool for Code Generation and Aspects Weaving[END_REF] code generation, which uses mapping rules scripts to produce source code files for a given target platform from UML models annotated with the MARTE profile 4 . Therefore, GenERTiCA is capable of generating code for many distinct languages (Java, C/C++, etc.), since there are mapping rules for the target platforms. The process is the same to generate code for different languages and therefore it is considered generic. This work proposes an extension for the GenERTiCA tool in terms of a new set of mapping rules to map UML meta-model elements into VHDL constructs. Following, a brief description of each step of AMoDE-RT design flow is presented.

Requirements Analysis and Identification

In the first phase, requirements and constraints of the distributed embedded real-time system are gathered. To accomplish this, the RT-FRIDA [START_REF] Freitas | DERAF: A high-level Aspects Framework for Distributed Embedded Real-Time Systems Design[END_REF] requirements analysis is performed, resulting in a set of documents describing system's requirements, functionalities and constraints. Afterwards, use case diagrams are created, depicting all expected functionalities of the distributed embedded real time system, and also the external elements that interact with the system.

Modeling

The next step is to specify the elements to handle the functional and non-functional requirements gathered in the previous phase. To model functional requirements, designers use class diagrams to describe the structure, and sequence diagrams to describe the methods behavior. Other behavioral and structural diagrams, such as activity or state diagrams, or composite structure or deployment diagrams, can also be used. However the class and sequence diagrams are mandatory to describe the structure and behavior of all system with correctness. These UML diagrams are annotated with the stereotype of the MARTE profile to specify real-time characteristics of (some) DES' elements. During this phase, the non-functional requirements handling are specified using aspects from the Distributed Embedded Real-Time Aspects Framework (DERAF) [START_REF] Wehrmeister | GenERTiCA: A Tool for Code Generation and Aspects Weaving[END_REF]. These aspects are modeled in the Aspects Crosscutting Overview Diagram (ACOD) [START_REF] Wehrmeister | An Aspect-Oriented Approach for Dealing with Non-Functional Requirements in Model-Driven Development of Distributed Embedded Real-Time Systems[END_REF], and the points (in the UML model) in which DERAF aspects perform adaptations are specified using Join Point Designation Diagrams (JPDD).

UML-to-DERCS Transformation

At this point GenERTiCA transforms the system specification, i.e. the UML model, into another model called DERCS (Distributed Embedded Compact Specification) [START_REF] Wehrmeister | GenERTiCA: A Tool for Code Generation and Aspects Weaving[END_REF], which represents an embedded system PIM free of information overlapping 5 . A UML specification can contain several model elements, representing the same element that hinders the code generation process. Thus these ambiguous elements of UML model are mapped in a single DERCS element, eliminating such ambiguities that could result in code with errors. When an inconsistency is detected, the UML-to-DERCS transformation algorithm stops and GenERTiCA informs this occurrence to the designer, requesting his/her intervention to solve the issue. Interested readers are referred to [START_REF] Wehrmeister | GenERTiCA: A Tool for Code Generation and Aspects Weaving[END_REF] to have more details on this UML-to-DERCS model transformation.

Code Generation

In this phase the code generation process executes a set of scripts (i.e. mapping rules), which guide the GenERTiCA tool to perform the model-to-text transformation from DERCS elements to constructions in the target platform. Furthermore, the code generation process also performs the aspects weaving. If the element under evaluation is affected by an aspect, the aspects weaving process modifies the generated code fragments according to aspects adaptations described in the mapping rules.

Mapping UML to VHDL

As previously mentioned, to generate code from the UML model, GenERTiCA adopts a script-based approach, in which small scripts define how to map model elements into target platform constructions, generating source code fragments that are merged to produce source code files. The script-based code generation process improves separation of concerns in mapping rules specification, as each script is concerned with the transformation of a single model element (or few of them) into source code fragment. In this sense, this work has proposed a set of mapping rules to allow VHDL code generation from UML models, following GenERTiCA's approach. Table 1 shows the mapping from UML concepts into VHDL ones. Scripts to accomplish these transformations have been developed and inserted in an eXtensible Markup Language (XML) 6 file, which guide GenERTiCA in the code generation process. Details on the created mapping rules are provided in the following sub-sections.

Classes are mapped into VHDL entity-architecture pairs. The class parameters are mapped to VHDL generic statements, while public attributes to VHDL entity ports and private attributes to VHDL signals. The methods are mapped to VHDL processes. The composition relationship, which describes the composition of a system from components, is mapped to a VHDL port map statements. Objects are instantiated as component structures into other entities. Events and Messages exchanges are implemented as entities ports that allow the communication between different entities and their processes (methods). Associations between classes are similar to the approach used to messages exchange, however component structures representing each associated-class (i.e. entities) are instantiated into the pair-class to accomplish the association by the mapping of signals between these two classes. Inheritance is obtained making the entity or architecture declaration with the VHDL key word "tagged", which means that the declaration is valid but still incomplete. Then, we declare a new entity/architecture using the VHDL key word "new", with the desired modifications. Static polymorphism is obtained using the VHDL configuration structure to bind the same VHDL component to different entities or architectures.

Mapping Rules

Mapping rules are specified as small scripts that create source code fragments (representing target platform constructions) from DERCS model elements. Source code files are made up from these generated code fragments. Scripts are stored and organized in one mapping rules file specified using the XML format. This XML file has a portable format allowing the specification of self-described content organized in a tree structure. XML tree organization facilitates scripts storage in terms of platform mapping rules repositories. It allows scripts to be reused in further projects that use the same target platform. Hence, the design effort to derive system implementation from an UML model is decreased. Scripts are located in the leaf nodes of the tree. The correct script is selected based on which element is being accessed by the code generation algorithm (i.e. the leaf node must match with the DERCS element). The code generation algorithm first identifies the type of the element. Afterwards, it tries to find the tree's leaf that better represents the type of such element being evaluated. Then, it executes the script contained within the found node (i.e. leaf), which will create a VHDL code fragment to that element being evaluated.

The language used to describe the scripts is the well-known open source scripting framework called Velocity7 , which defines the Velocity Template Language (VTL) that provides all functionalities required to assist the code generation approach implementation. VTL is a Java-like scripting language, which returns a string as result of script execution. Thus, the generated source code fragment is obtained by means of accessing model information through DERCS API.

To illustrate what is a script an example is given above. This script is responsible for the code generation to the classes' methods. It generates one VHDL process for each method in the classes.

Concepts

The concepts used in this work to the transformation of UML structures into VHDL code are based in [START_REF] Mcumber | UML-Based Analysis of Embedded Systems Using a Mapping to VHDL[END_REF] and [START_REF] Ecker | An object-oriented view of structural VHDL description[END_REF]. These works gave us some ideas on how to represent in VHDL the structure of classes, attributes, methods, association between classes, events and messages exchange, inheritance, and polymorphism. UML Model is object-oriented, while the VHDL code is structured. This semantic gap between abstraction levels hinders the mapping between these two languages. It may be one of the main reasons why it is very rare, until the present days, to find works addressing to VHDL code generation from UML models and/or any commercial tool.

From these concepts some rules have been developed and tested to this version of the VHDL mapping rules. They represent the first version of the mapping rules to provide the VHDL code generation through GenERTiCA tool. Up to now, it has been developed and tested mapping rules to generate VHDL code from UML classes, attributes and behaviors. Their feasibility is shown in the case study of next section. The concepts of inheritance, polymorphism, associations between classes, etc, are being be implemented in the next version of this work.

Case Study

This section shows an example of automatic VHDL code generation from a UML model. The system under evaluation is composed by an automatic valve and the sensors that give information about the valve's states. This valve is used to regulate the water flow and is part of CISPI (Conduite de grands systèmes industriels à risque) experimental platform located at CRAN (Centre de Recherche en Automatique de Nancy -France). By applying the proposed approach, we intend to integrate new functionalities supported by a FPGA in this valve, leading to an implementation of an intelligent component. This intelligent component is part of a mechatronic system which also contains other mechanical, electronics and computational parts. The electronics and computational parts represent the control system, which is composed of logical functions executing in a hardware platform. Logical functions represent components' behavior and are usually implemented in software. This work aims at implementing these logical functions as hardware. Thus the system has been specified in UML and its implementation has been generated as VHDL code. The generated code represent the hardware description of the logical and the control system parts (logical functions), which is executed inside of a FPGA (hardware platform).

The development of this case study was compliant with the AMoDE-RT approach defined in the section III. Then, the first step was to gather the requirements and create a use case diagram to identify system's services and the actors that interact with the system. Only two actors were considered. The "User" actor request services to open and close the valve, while the "Maintenance Operator" actor requests services to know the number of times the valve performed opening and closing actions. Information about these numbers is used for assessing the component physical heath.

The second step was to create class and sequence diagrams. The class diagram has been built from the knowledge acquired in the use case. It represents the system's structure. All services have been modeled as classes. These classes work together to provide the system's services. They are enclosed by a main class responsible for all system. No DERAF aspects have been used to deal with non-functional requirements, and hence, no ACOD and JPDD diagrams have been created.

Services are modeled as sequences of actions in sequence diagrams. Sequence diagrams represent the exchange of messages between the objects that compose the system in order to represent the expected system behavior. Each service demanded by the actors results in the execution of one or many sequence diagrams. Fig. 2 depicts a part of the main sequence diagram of valve's system related to the service solicitations. The rule implementing the association between classes was not developed in this first version of the mapping rules. Thus, all this case study was re-modeled as just one class and one method. From these new diagrams the feasibility of the system could be tested and the 3 rd and 4 th steps performed. The mapping has transformed the UML class and sequence diagrams into one single VHDL file, which contains an entityarchitecture pair declaration. The UML class is represented by the entity-architecture pair, while the UML behavior (method) is represented by a VHDL process inside of the architecture.

The code that appears in fig. 3 is the result of the automatic code generation performed by the mapping rules proposed in this work. This first test generated code for only one embedded component. For this component have been automatically generated 80 lines of VHDL code, covering 100% of the needed code for the application. This component is simple and its model was developed on only one UML class, however the results are encouraging for future works.

The real feasibility of the VHDL code can be tested quickly, since this code can be synthesized in the resulting bitstream uploaded in a FPGA development platform, such as Virtex-II PRO (V2-Pro) development system by Digilent 8 .

Related Works

This section discusses some projects and commercial tools that propose transformations from UML specifications to source code VHDL. Among these works different approaches to generate source code from UML models have been found. Some use only one diagram (e.g. class diagram) to generate code, while others use a combination of distinct diagrams (e.g. class and state diagrams, sequence and/or activities diagrams) [START_REF] Long | Consistent Code Generation from UML Models[END_REF]. Thus the presented related works generate code ranging from classes skeletons to code containing system elements behavior.

In [START_REF] Mcumber | UML-Based Analysis of Embedded Systems Using a Mapping to VHDL[END_REF], a framework has been developed to derive VHDL specifications from UML's class and state diagrams. They use homomorphic mappings between dissimilar structures while preserving metamodel class associations in a way that resembles MDA technique. However, their generated VHDL code focuses on simulation and verification of UML models rather than on hardware synthesis.

An interesting work has been developed in [START_REF] Coyle | From UML to HDL: a Model-Driven Architectural approach to hardware-software co-design[END_REF], in which the MODCO tool is presented. It uses MDA techniques to define high-level model-based system descriptions that can be implemented in either hardware or software. Thus it can transform UML state diagrams directly into synthesizable VHDL. State machines in UML are used to describe the behavior of a part of a system. However, the complete code is generated for the behavior. Functional requirements are mapped to UML component, class, use case and state diagrams. Non-functional requirements are specified as UML annotations that describe performance constraints using propertyvalue pairs defined by UML profiles. However, the authors targeted flat statetransition diagrams without supporting hierarchy and concurrency, and also only covering a small subset of UML state diagram constructs.

In [START_REF] Wood | A Model-Driven development Approach to Mapping UML State Diagrams to Synthesizable VHDL[END_REF], they have developed a framework for deriving VHDL specifications from UML state diagrams, and also a set of rules, which enable automated generation of synthesizable VHDL code from UML. Their engineering process is based on metamodels. Concepts of the UML state diagram metamodel are mapped onto concepts of the VHDL metamodel. There are two transformations between models, happening in the following way: the first transformation converts the main UML model into state diagram models; and the second one maps the state diagram models onto concepts in the VHDL language. A model-to-text transformation is used to generate synthesizable VHDL code from the VHDL model.

Two commercial tools to generate VHDL code could also be found. StateCAD9 by Xilinx is a graphical entry tool for digital design that has its own graphical notation to represent state diagrams as bubble diagrams. StateCAD automatically generates HDL (VHDL and Verilog) code, for simulation and also synthesis, directly from these state diagrams. The other tool is Simulink HDL Coder10 by MathWorks, which generates synthesizable VHDL and Verilog code from simulink models, stateflow charts, and embedded Matlab code. Simulink HDL Coder also generates simulation and synthesis scripts, enabling to simulate and synthesize quickly the developed design. These commercial tools do not generate VHDL code from UML specifications.

The approach proposed in [START_REF] Wood | A Model-Driven development Approach to Mapping UML State Diagrams to Synthesizable VHDL[END_REF] realizes the mapping between models, similarly to other MDE techniques. In [START_REF] Coyle | From UML to HDL: a Model-Driven Architectural approach to hardware-software co-design[END_REF], there is a separation of the functional and nonfunctional requirements in the modeling stage using distinct UML diagrams. These mentioned works are limited, because they cover a specific subset of the UML structures, and also use only the UML state diagrams. By using GenERTiCA the designer can use distinct UML diagrams, combining them to specify the full functionality in terms of structure, behavior and non-functional requirements handling, since some guidelines are followed. These guidelines are simple and intuitive, allowing designers to separate the functional and non-functional requirements. Moreover, by using the extension proposed in this work, it is possible to generate code for VHDL and also for the Java and C++ languages from the same model.

Conclusions

This work addresses the problem of generating HDL descriptions from UML models. It presented the initial set of mapping rules to generate VHDL code from class and sequence diagrams, using GenERTiCA tool. To achieve this goal, is has been proposed a mapping from object-oriented concepts supported in UML into concepts used by VHDL. Then, a set of mapping rules (used by GenERTiCA to generate VHDL code) has been developed. These rules extend the functionality of GenERTiCA tool, allowing it to generate HDL code in addition to software source code. Besides, this paper has described all steps of the proposed approach that must be followed to generate automatically VHDL descriptions.

To demonstrate the proposed approach, a case study has been presented. It showed a small part of a maintenance system, i.e. automatic control of a valve implemented as a Smart Component. As mentioned mapping rules have been implemented and tested to produce VHDL code from UML's classes, attributes and behavior. Results shown that, for developing simple systems, 100% of the necessary code could be generated. Hence, despite the case study's size, the results are considered satisfactory since we see great potential to scale the approach to more complex systems.

To continue this work the following direction will be pursued: to complete the rules needed for the code generation of VHDL structures. The concepts of inheritance, polymorphism, associations between classes are very important to be able the modeling of complex systems; to develop the rules for the non-functional requirements implemented by aspects; to test and to prove the new rules; to perform more tests with FPGA boards; to apply the proposed approach in a complex real system, such as an industrial maintenance and prognostic systems.

Fig. 1 .

 1 Fig. 1. General overview on the whole process.

Fig. 2 .

 2 Fig. 2. Sequence diagram of the main function in the valve's system.

Fig. 3 .

 3 Fig. 3. VHDL code representing the architecture declaration and run method initial part.

Table 1 .

 1 Concepts Mapping.

	UML Element	VHDL Element
	Class	Entity-Architecture pair
	Public attribute	Entity Ports
	Private attribute	Signals
	Methods	Processes
	Events and Message exchange	Entity Ports
	Associations between classes	Entity Ports
	Inheritance	VHDL key word "new"
	Static polymorphism	Configuration structure
	Objects instantiation	Component structure

UML

2.2 Specification. Object Management Group, http://www.omg.org/ spec/UML/2.2/ 2 http://www.ibm.com/software/rational/ 3 http://www.artisansoftwaretools.com

UML Profile for Modeling and Analysis of Real-time and Embedded Systems (MARTE), http://www.omg.org/spec/MARTE/1.0

Information overlapping in UML models means the same feature of the target system, which has been specified using distinct diagrams depicting different viewpoints of the same structural/behavioral characteristic.

eXtensible Markup Language (XML) 1.0 (Fifth Edition). http://www.w3.org/TR/2008/REC-xml-20081126/

http://velocity.apache.org

Digilent Inc. www.digilentinc.com

http://www.xilinx.com/

http://www.mathworks.com/products/slhdlcoder/


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