Héctor Fernández
email: hector.fernandez@inria.fr

Cédric Tedeschi
email: cedric.tedeschi@inria.fr

Thierry Priol
email: thierry.priol@inria.fr

A Chemistry-Inspired Workflow Management System for Scientific Applications in Clouds

Keywords: Scientific workflows, Workflow management system; Nature-inspired models, Chemical programming

With the proliferation of Web Services, scientific applications are more and more designed as temporal composition of services, commonly referred to as, workflows. To address this paradigm shift, different workflow management systems have been proposed. If their efficiency has been established over centralized reliable systems, it is questionable over highly decentralized failure-prone platforms. Scientific applications started to be deployed over emerging clouds, leading to new issues, like elasticity, i.e., the possibility to dynamically refine at runtime the amount of resources dedicated to an application. This raised a new demand for programming models, able to express autonomic self-coordination of services in a dynamic, elastic platform.

Chemistry-inspired computing recently regained momentum in this context, naturally expressing parallelism, distribution, and autonomic behaviors. While its high expressiveness and adequacy for this context has been established, the chemical model severely suffers from a lack of proof of concepts. In this paper, we concretely show how to leverage such models in this context. We focus on the design, the implementation and the experimental validation of a chemistry-inspired scientific workflow management system.

I. INTRODUCTION

Until recently, scientific applications were commonly hardto-maintain unreadable scripts, leading to a poor reusability level and high maintenance costs. With the proliferation of Web services and the increasing adoption of service-oriented computing, whose primary goal is to make a collection of software services accessible through the network, scientists started to develop their applications as compositions of Web services, today commonly referred to as workflows. This shift of paradigm recently led to more reutilization, and experiment sharing in the community. The specification and execution of such workflows are managed by workflow management systems, responsible for the coordination of the services involved. Addressing the limitations of initial workflow languages such as the BPEL business standard [START_REF] Oasis | Web services business process execution language, (WS-BPEL)[END_REF], different systems, for example Taverna [START_REF] Oinn | Taverna: lessons in creating a workflow environment for the life sciences: Research articles[END_REF], Kepler [START_REF] Ludascher | Scientific workflow management and the kepler system: Research articles[END_REF], Triana [START_REF] Taylor | The Triana Workflow Environment: Architecture and Applications[END_REF], Pegasus [START_REF] Deelman | Pegasus: A framework for mapping complex scientific workflows onto distributed systems[END_REF] or Askalon [START_REF] Fahringer | Askalon: a tool set for cluster and grid computing: Research articles[END_REF] provide nice features such as implicit parallelism or data-driven coordination, increasing the abstraction regarding execution management while improving the efficiency and manageability of science workflows, as formulated by Zhao, Raicu, and Foster in 2008 [START_REF] Zhao | Scientific workflow systems for 21st century, new bottle or new wine[END_REF].

With the rise of cloud computing, science workflows started to be deployed over clouds. For example, the Magellan project is aimed at providing cloud infrastructures for science [START_REF]The magellan research project[END_REF]. Clouds, or not far from here, federation of clouds, are the new target platform for scientific workflows. One key properties of cloud computing is elasticity, i.e., the possibility for cloud users, to dynamically adapt the quantity of resources at his disposal at runtime.

Thus it appears that future scientific workflow systems and languages should provide a natural way to express both workflows and platform characteristics. We identify several critical features the future WMS must address: (i) the high degree of parallelism and distribution of services deployed over a federation of clouds, (ii) the potential issues brought by a centralized coordinator such as single point of failure and scalability, and (iii) the elasticity and distribution of clouds, especially when dealing with federations of clouds.

Lately, nature metaphors, and in particular chemistryinspired analogies have been identified as a promising source of inspiration for developing new approaches for autonomous service coordination [START_REF] Viroli | A biochemical approach to adaptive service ecosystems[END_REF]. Among them, the chemical programming paradigm is a high-level execution model. Within such a model, a computation is basically seen as a set of reactions consuming some molecules of data interacting freely within a chemical solution producing new ones (resulting data). Reactions take place in an implicitly parallel, autonomous, and decentralized manner. More recently, the Higher-Order Chemical Language (HOCL) [START_REF] Banâtre | Generalised multisets for chemical programming[END_REF] raised the chemical model to the higher-order, providing a highly-expressive paradigm: every entity in the system (in our case data, services and their dependencies, and the platform itself) is seen as a molecule. Moreover, rules can apply on other reaction rules, programs dynamically modifying programs, opening doors to dynamic adaptation. This model is now envisioned as an alternative naturally expressing autonomous coordination [START_REF] Banâtre | Chemical Programming of Future Service-oriented Architectures[END_REF]. However, while its expressiveness and adequacy to service coordination have been established, the actual experimentation of the chemical model has remained quite limited until now. There is a strong need of a proof of concept to show its viability, in particular compared to current reference WMS.

Contribution. In this paper, we present a workflow managment system able to solve a wide variety of workflow patterns both in a centralized and a decentralized way following the chemical model. Its implementation and performance evaluation on different classic scientific workflows presenting different characteristics are discussed. For the sake of comparison and discussion, workflows tested were also executed on top of Taverna and Kepler WMS, validating our software prototype, and establishing the viability of the concept.

Section II introduces the preliminaries of our work: workflow management systems and the chemical computing paradigm. Section III describes the architecture and workflow engine we have built. We show how a workflow is described in our model, and how chemical rules are defined and combined so that they can solve a wide variety of workflow patterns. Section IV focus on the implementations of both centralized and decentralized versions of the system defined. Section V details the experimental campaign and its results. Section VI discusses related works. Section VII draws a conclusion.

II. BACKGROUND

We now introduce the two main areas of this work: workflow management systems and the chemical computing model.

A. Workflow Management for e-Science

The increase in the reliance on service-oriented architectures (SOA) in e-science resulted in applications to be more and more defined as workflows of services. As a natural consequence, workflow management systems have gained recently considerable attention. The BPEL standard [START_REF] Oasis | Web services business process execution language, (WS-BPEL)[END_REF] and followers [START_REF] Barker | Scientific Workflow: A Survey and Research Directions[END_REF] were first briefly adopted by the scientific community. Then, science-oriented workflow languages and systems were designed, to cope with the different characteristics of scientific applications, in parallelism high parallelism and need for scheduling. In this way, a number of systems were designed for the expression and execution of scientific workflows. Taverna [START_REF] Oinn | Taverna: lessons in creating a workflow environment for the life sciences: Research articles[END_REF], Kepler [START_REF] Ludascher | Scientific workflow management and the kepler system: Research articles[END_REF], Triana [START_REF] Taylor | The Triana Workflow Environment: Architecture and Applications[END_REF], Pegasus [START_REF] Deelman | Pegasus: A framework for mapping complex scientific workflows onto distributed systems[END_REF] or Askalon [START_REF] Fahringer | Askalon: a tool set for cluster and grid computing: Research articles[END_REF] provide nice features such as implicit parallelism or data-driven coordination, increasing the abstraction regarding execution management while improving the efficiency and manageability. All typically provide a visual notation for composition of services.

The remainder of this section introduces the two opensource workflow systems used for the sake of validation of our work: Taverna and Kepler. We chose Kepler and Taverna because they are among the most used and mature open-source scientific systems. As mature as it can be, we choose not to use Pegasus as the resource management is integrated into the workflow manager, which is not of our primary concern here.

Taverna builds upon service-oriented architectures, and the Web service standards. Interactions between services (processors) are defined using the XML-based Scufl language or a GUI. Taverna is data-driven, data-dependencies specify links among different services, thus parallelism is implicit. Note that control-dependencies can also be specified by links that define precedence conditions among processors. Taverna's workflow engine is centralized; a unique coordinator manages the coordination of all computation blocks.

Kepler is a centralized workflow engine built upon Ptolemy II [START_REF] Buck | Ptolemy: A framework for simulating and prototyping heterogenous systems[END_REF] by the members of the Science Environment for Ecological Knowledge (SEEK) project and the Scientific Data Management (SDM) project. It also relies on a data-driven model for simulating, and designing of real-time and concurrent workflows using a proprietary modeling markup language called MoML. This language is based on the actor-oriented modeling paradigm which consists in a composition of computation blocks called actors that represent operations or data sources. Due to its data-driven behavior, Kepler provides an intuitive and implicit parallel execution, however it hinders the execution of complex workflow patterns. In recent versions, some control structures can be supported through somewhat sophisticated programming.

The limitation of both Taverna and Kepler is the lack of (i) facilities to describe more complex control-flow patterns, and (ii) support for decentralized coordination of the workflow execution.

B. Chemical Programming

Nature analogies, and more specifically bio-chemical metaphors, have recently gained momentum in the construction of programming models coping with the requirements of the Internet of Services [START_REF] Viroli | A biochemical approach to adaptive service ecosystems[END_REF]. Initially proposed to naturally express highly parallel programs, the chemical programming paradigm exhibits properties required in emerging service platforms and naturally express autonomic coordination.

According to the chemical metaphor, molecules (data) float in a chemical solution, and react according to reaction rules (program) producing new molecules (resulting data). These reactions take place in an implicitly parallel, autonomous, and non-deterministic way until no more reactions are possible, a state referred to as inertia. The computation is carried out according to local conditions without any central coordination, ordering or serialization. This programming style allows writing programs cleared of any artificial sequentiality and to concentrate on the functional aspects of the problem solved. The presence of a molecule is enough to trigger a reaction requiring such a molecule. Nevertheless, as it will be shown, we can express sequentiality if needed.

The GAMMA language [START_REF] Banâtre | The gamma model and its discipline of programming[END_REF] first formalized chemical programming. In GAMMA, the solution is a multiset containing the molecules, and reactions between molecules are rules rewriting the multiset. The multiset is the unique data structure. More recently, higher-order chemical programming has been proposed, through HOCL (Higher Order Chemical Language) [START_REF] Banâtre | Generalised multisets for chemical programming[END_REF]. In HOCL, every entity is a molecule, including reaction rules. A program is a solution of molecules, formally a multiset of atoms, denoted A 1 , A 2 , . . . ,A n , "," being the associative and commutative operator of construction of compound molecules. Atoms can be constants (integers, booleans, etc.), reaction rules, tuples of n atoms, denoted A 1 :A 2 :. . .:A n , or sub-solutions, denoted M i , where M i is the molecule content of the sub-solution. A reaction involves a reaction rule replace P by M if V and a molecule N satisfying the pattern P and the reaction condition V . The reaction consumes the molecule N to produce a new molecule M. This rule can react as long as a molecule satisfying the pattern P exists in the solution. Its one-shot variant, denoted one P by M if V , reacts only once, and is consumed in the reaction. Rules can be named or appear explicitly in the solution. Let us consider the simple HOCL program below that extract the maximum even number from a set of integers. The selectEvens rule removes odd numbers from the solution, by repeated reactions with an integer x, ω denoting the whole solution in which selectEvens floats, deprived of x. The getM ax rule reacts with two integers x and y such that x ≥ y and replaces them by x. In a solution of integers, this rule, by its repeated application extracts the maximum value. The solution is composed by (i) a sub-solution containing the input integers along with the selectEvens rule, and (ii) a higher-order rule (on Line 1.06) that will open the subsolution, extract the remaining (even) numbers and introduce the getM ax rule.

Solving the problem requires the sequentiality of the reactions of the two rules. This can be achieved by the higherorder: in an HOCL program, a sub-solution can react with other elements only if it has reached the state of inertia. In other terms, the higher-order rule will react with the subsolutions only when no more reactions are possible within it, i.e., when it contains only even numbers. (Note that the order in which odd numbers are deleted is non-deterministic.) The result is then as follows:

selectEvens, 2, 6, 8 , replace-one selectEvens = s, ω by getM ax, ω Then, the higher-order rule reacts with it, extracting remaining numbers, introducing dynamically the getM ax rule, and in this way triggering the second phase of the program where the maximum value is kept. The resulting solution is: 2, 6, 8, getM ax . getM ax then reacts with pairs of integers until only 8 remains. Note that, due to the higher-order, putting both rules directly in the solution of integers could entail a wrong behavior as the pipeline between the two rules would be broken, possibly leading to a wrong result, for instance if getM ax reacts first with the 8 and 9, thus deleting the 8.

While this example is quite simple, it already provides the intuition behind autonomic coordination and online adaptation.

III. CHEMICAL WORKFLOW MANAGEMENT SYSTEM

In this section, we describe chemistry-inspired workflow management systems. First, the coordination mechanisms developed, which build upon higher-order chemistry, are presented. Then, the architecture underlying it, for both centralized or decentralized coordination are described. The concepts presented in this section take their origin in the founding work presented in [START_REF] Fernández | Decentralized Approach for Execution of Composite Web Services using the Chemical Paradigm[END_REF].

A. Workflow Representation

We now consider a workflow of Web Services (WS). The general shape of the chemical representation of a workflow is as follows: the main solution is composed of as many subsolutions as there are WSes in the workflow. Each sub-solution represents a WS with its own data and control dependencies with other WSes. More formally, a WS is a molecule of the form W Si : . . . where WS i refers to the symbolic name given to the service whose connection details and physical position are hidden.

Let us consider a simple workflow example whose chemical representation is illustrated by Figure 1. It is composed of four services S 1 , S 2 , S 3 and S 4 . In this example, after S 1 completes, S 2 and S 3 can be invoked in parallel. Once S 2 and S 3 have both completed, S 4 can be invoked. W S1 : . . . to W S4 : . . . represent WSes in the solution. The relations between WSes are expressed through molecules of the form DEST:WSi with W Si being the destination WS where some information needs to be transferred. For instance, we can see in the W S1 sub-solution that WS1 will transfer some information (its outcome) to WS2 and WS3 (Line 2.02). Let us have a more precise look on these dependencies. WS2 has a data dependency: it requires a molecule RE-SULT:WS1:value1 containing the result of S 1 to be invoked (second part of Line 2.04). The two molecules produced by the reaction represent the call to S 2 and their input parameters. They are expressed using a molecule of the form CALL:Si, and a molecule PARAM: in 1 ,...,in n , where in 1 ,...,in n represent the input parameters to call the service Si. In Figure 1, this input parameter corresponds to the result of some previous service S j . WS3 works similarly. WS4 performs a particular control pattern known as synchronization. It needs to wait until WS2 and WS3 have completed, in other words, until the molecules RESULT:WS2:value2 and RESULT:WS3:value3 appear in its sub-solution, to start its own execution. In addition, a data dependency is also expressed in WS4: the result of S 2 is required to call S 4 .

To ensure the execution of a chemical workflow, additional generic chemical rules (i.e., independent of any specific workflow) must be defined. These rules consume and generate additional molecules to manage transfer of information between services, condition checking, fault detection, and more complex control flows. To express the whole logic of a workflow, these rules have to be composed relying on the analogy of molecular composition. This concept consists in the composition of several molecules, which are composed based on data molecule dependencies, and whose reactions produce new molecules reacting in their turn, and so on, to complete the workflow.

B. Generic Rules for Invocation and Transfer

Common tasks in a workflow of services are service invocation and information transfer between services. We now review three generic rules illustrated in Algorithm 1, responsible for these tasks, and that will be commonly encountered in the compositions presented later. The invokeServ rule encapsulates the actual invocation of services. When reacting, it invokes the Web Service S i , by consuming the tuples CALL:Si representing the invocation itself, and PARAM: in 1 ,...,in n representing its input parameters, and generates the molecules containing the results of the invocation in the WSi sub-solution. The molecule FLAG INVOKE is a flag whose presence in the solution indicates that the invocation can take place. The preparePass rule is used for preparing the messages to transfer the results to their destination services, that will later trigger the execution of the passInfo rule.

Algorithm 1 Basic generic rules. Rule passInfo transfers molecules of information between WSes. This rule reacts with a molecule WSi: PASS:d: ω 1 that indicates that some molecules (here denoted ω 1 ) from WSi needs to be transfer to d. These molecules, once inside the sub-solution of d will trigger the next step of the execution. Therefore, the molecule ω 1 will be transferred from subsolution WSi to sub-solution d, when reacting with passInfo rule.

C. Complex Workflow Patterns

With generic rules described until now, the engine can only support data flows. However, more complex control flows are required to be taken into account, in order to solve a broader range of workflow definitions. We now illustrate how HOCL can be leveraged to deal with complex control flows, by detailing a particular pattern known as Simple Merge.

As illustrated in Figure 2, a simple merge pattern is similar to a XOR operation. It involves a structure where two or more source service flows converge into a single destination in an asynchronous way. The destination service must however be launched only once, regardless of the number of incoming. In other words, only the first source service to complete will influence the remainder of the workflow execution. To enhance our workflow engine with the support of the simple merge pattern, we need to define the appropriate generic rules and dispatch them in the sub-solutions of WSs involved. These rules are given in Algorithm 2. The sm preparePass reaction rule is used to add, in the subsolution of every incoming service, a particular MERGE:Yes molecule to the information to be transferred to the destination service (see Lines 4.01 and 4.03). The destination WS waits for this molecule and only the first MERGE:Yes molecule received in its sub-solution will be consumed. Next, sm setFlag reaction rule takes place producing one molecule of the form FLAG INVOKE, allowing the service invocation. The following MERGE:Yes molecules received will be ignored. In terms of molecular composition, each source WS will have in its sub-solution one sm preparePass rule (A on Figure 2) and one passInfo (B on Figure 2) rules, they are composed with sm setFlag rule (C) in the destination WS. For space reasons, we omit more complex control flows (such as synchronization merge, exclusive choice or discriminator). However, the support for a wide range of control flow patterns can be found in the research report [START_REF] Fernández | Self-coordination of Workflow Execution Through Molecular Composition[END_REF], as well as its design process.

To sum up, the coordination is achieved through a set of autonomic, and local reactions taking place within each WS's sub-solution (or between two WSs' sub-solutions), providing all the abstractions for a natural expression of a decentralized coordination of virtually all identified workflow patterns [START_REF]The workflow patterns website[END_REF].

D. Architectural Frameworks

We now show how the chemical engine can be powered over both centralized and decentralized architectures.

1) Centralized Architecture: Following the examples of most of Workflow Management Systems mentioned before, the coordination can be managed by a single node, referred to as the chemical workflow service, as illustrated by Figure 3. First, notice the S components. They represent the interface with the actual distant WSes to be called. Then, the workflow definition and coordination information as described before (i.e., the multiset) is accessed by the chemical engine that will perform the reactions required. 2) Decentralized Architecture: Distribute the control means that each service involved will take its part in the coordination process. Each WS is now chemically encapsulated, to form a Chemical Web Service (ChWS). Each ChWS is equipped with a chemical engine and a local copy of part of the multiset on which its chemical interpreter will act, to realize its part of the coordination. The multiset, containing the workflow definition and thus all required coordination information, will now act as a space shared by all ChWSes involved in the workflow. In other words, ChWSes will communicate by reading and writing it, as illustrated by Figure 4.

IV. SOFTWARE PROTOTYPE

In this section, we discuss the implementation of a software prototype for both previously described architectures. The low layer of our prototype is an HOCL interpreter based on on-thefly compilation of HOCL specifications. The whole prototype is written with Java.

A. Centralized Version

The prototype is illustrated by Figure 5. As mentioned in Section III, the workflow definition is executed as a chemical program by the chemical workflow service. The low layer of the architecture is an HOCL interpreter. Given a workflow specification as input (an HOCL program), it executes the workflow coordination by reading and writing the multiset initially fed with the workflow definition. The interface between the chemical engine and the distant services themselves is realized through the service caller. The service caller relies on the DAIOS framework [START_REF] Leitner | Daios: Efficient dynamic web service invocation[END_REF], which provides an abstraction layer allowing dynamic connection to different flavours of services (SOAP or RESTFul), at the same time abstracting the target service's internals. DAIOS was specially extended with a module which automatically generates the bindings, as well as input and output messages required between the chemical engine and a Web service.

B. Decentralized Version

The decentralized prototype is illustrated in Figure 6. Basically, the difference between centralized and decentralized implementations is the functionality of the multiset. It now represents a shared space playing the role of a communication mechanism and a storage system. As we detailed before, the workflow definition is comprised of one sub-solution per WS involved. The information in one sub-solution can only be accessed by the ChWS owner of/represented by that sub-solution. On each ChWS, a local storage space acts as a temporary container for the sub-solution to be processed by the local HOCL interpreter. ChWSes communicate with the multiset through RMI interfaces to read and modify it concurrently. Periodically and independently from each other, ChWSes request their sub-solution to this multiset. The subsolution obtained is then processed, modified by the local HOCL interpreter and then pushed back to the multiset as an update. On the other hand, the service caller now represents the interface between a ChWS and a concrete WS. To implement the needed communication mechanism, some Java RMI modules are included into the ChWSes and the multiset. The multiset has been enhanced with an RMI servercommunicator module, providing three remote operations for the ChWSes: put operation adds molecules in the sub-solution of a specified ChWS, get withdraws the content of a ChWS solution, and isConsumable determines whether there is any new information to be processed by the calling ChWS. On the ChWS' side, an RMI listener module has been integrated on each ChWS, periodically calling the isConsumable operation, which, if it returns true, triggers the invocation of the get operation to extract the content of the WSi solution which will be copied into its local multiset. Once the HOCL interpreter has finished the execution, internally the RMI listener call the put operation to update the content of its sub-solution into the multiset. This multiset implementation presents some similarities with tuplespaces, founded in the Linda coordination language [START_REF] University | Generative communication in linda[END_REF].

V. PERFORMANCE RESULTS

In this section, we present and analyse our experimental results. Four engines have been used: Taverna Workbench 2.2.0, Kepler 2.0, and both centralized and decentralized version of the HOCL-based workflow engine, referred to as HOCL Cen. and HOCL Dec. in the following, respectively.

Recall that our objective is not so much comparing performance as is establishing the viability of a chemistry-based workflow engine. In other terms, Taverna and Kepler represent validated standards we use as guidelines.

A. Workflows considered

Three scientific workflows were executed. Illustrated by Figure 7 (left), BlastReport3 is a bio-informatics workflow which retrieves a blast report of a protein in a database given its ID. The second one, CardiacAnalysis 1 , illustrated in Figure 7 (right), is a cardiovascular image analysis workflow which extracts the heart's anatomy from a series of image sequences by applying image processing algorithms, developed by the CREATIS-LRMN biomedical laboratory 1 . The third one, Montage [START_REF] Berriman | Montage: A grid enabled engine for delivering custom science-grade mosaics on demand[END_REF], given in Figure 8 is a classic astronomical image mosaic workflow processing large images of the sky. These workflows present different characteristics such as the number of services involved, the amount of data exchanged and the complexity of the coordination required (data processing included, such as iterations of lists of objects). We attempt to characterize these workflows as follows.

• The BlastReport workflow includes 5 services, and presents a medium level of data exchange (simple objects, lists) and low coordination overhead -it is composed mostly of sequences. • The CardiacAnalysis workflow includes 6 services, presenting a high level of data exchange (complex objects, lists) and a high coordination overhead (synchronizations, loop iteration, parallelism). This overhead does not appear on Figure 7 (right), and is due to the re-entrant nature of the services. Indeed, some services produce lists of The workflow definitions used for each WMS are available online 2,3 .

B. Centralized experiments

The workflows were first run using Taverna, Kepler, and HOCL Cen, on a local machine equipped with the Intel coreduo T9600 2.8 Ghz processor and 4GB of memory. Figures 9 and 10 present the results. A first encouraging result is that the execution time for the Montage workflow, i.e. for a workflow with limited data exchange and coordination overhead, on Kepler, Taverna and HOCL Cen. are quite similar, and even slightly reduced on the HOCL Cen. WMS.

For the BlastReport workflow, while results are again similar for the different WMSs, HOCL Cen. takes a little more time. This can be explained by the increased size of the multiset representing the BlastReport workflow (in terms of number of molecules). However, in terms of ratio, execution times remain very close. Finally, as we can see in Figure 10, the increased coordination overhead of the CardiacAnalysis. As mentioned before, this workflow relies on a lot of data processing related to the coordination itself, which, in the case of HOCL Cen., results in a significant increase of the size and processing time of the multiset. Also, no support for parallel execution has been implemented on our current HOCL interpreter version. These two optimization aspects will be investigated in the future.

C. Decentralized experiments

The workflow were also executed following the HOCL Dec. design. The experiments were conducted on the Grid'5000 platform [START_REF]Grid'5000[END_REF], specifically, on the paradent cluster, located in Rennes, each node being equipped with two quad-core Intel Xeon L5148 LV processors, 30 GB of RAM and 40GB InfiniBand Ethernet cards. We now focus on the right-most bar of each series of Figures 9 and10.

A first observation is the performance degradation using HOCL Dec. on the Montage workflow. While the coordination is executed locally on each ChWS (here the coordination is shared among 27 services), the multiset remains a space shared by every ChWSes leading to potential scalability issues. Also, RMI is known to have different some scalability weaknesses related to memory management [START_REF] Alt | Adapting java RMI for grid computing[END_REF], in particular provoked by the distributed garbage collector.

On the BlastReport, a performance gain over HOCL Cen. is obtained with HOCL Dec., thanks to the distribution of the coordination over the 5 services involved. The number of services participating is not high enough to provoke the previous scalability problems on the multiset side.

For the CardiacAnalysis workflow, a considerable performance gain is also obtained using HOCL Dec., demonstrating the benefits of a decentralized workflow execution when workflows present a high coordination overhead like Cardiac-Analysis. Exploiting the processing resources of each ChWS, the list handling and adaptation tasks are separately managed by each ChWS.

D. Discussion

This series of experiments leads to several conclusions. They constitute a proof of the viability of a chemistry-based engine, as for some workflows, its performance are similar to those of Kepler and Taverna.

Nevertheless, scalability issues need to be tackled. The first limitation comes from the design of the decentralized architecture itself, in which the multiset can constitue a bottleneck. To deal with the decentralization of the multiset itself, we recently formulated solutions based on peer-to-peer protocols, able to distribute and retrieve objects (here, molecules) at large scale [START_REF] Obrovac | On Distributing the Runtime of the Chemical Programming Model[END_REF]. One of the next steps of this work is to build the HOCL Dec. environment on top of such approaches to remove the bottleneck problem, and propose a fully decentralized workflow engine. A more short term possibility is to optimize the communication scheme used, currently based on RMI, also known for some scalability limitations. We will consider relying on a different communication tool.

Recall, beyond performance or optimization considerations, that the chemical models provide all the needed abstractions to naturally express both data-driven and complex controldriven execution, including particular features like cancellation. Please refer to [START_REF] Fernández | Self-coordination of Workflow Execution Through Molecular Composition[END_REF] for details. We consider the chemical abstraction as participating in the long term objective of improving the workflow execution models on emerging platform, like clouds, where the elasticity brings new modeling challenges.

VI. RELATED WORKS

Workflow management systems generated a high amount of literature and systems. With the proliferation of serviceoriented architectures, and the adoption of the Web service standard, several service-based WMS for scientific applications was proposed, such as Taverna [START_REF] Oinn | Taverna: lessons in creating a workflow environment for the life sciences: Research articles[END_REF], Kepler [START_REF] Ludascher | Scientific workflow management and the kepler system: Research articles[END_REF], Triana [START_REF] Taylor | The Triana Workflow Environment: Architecture and Applications[END_REF], Pegasus [START_REF] Deelman | Pegasus: A framework for mapping complex scientific workflows onto distributed systems[END_REF] or Askalon [START_REF] Fahringer | Askalon: a tool set for cluster and grid computing: Research articles[END_REF]. A common limitation of these systems is their lack of support for dynamic adaptation and decentralization, leading to scalability, as well as reliability issues. Our decentralized architecture share some similarities with the Linda [START_REF] University | Generative communication in linda[END_REF] approach, which pioneered the coordination through a shared space. In the same vein, recent works address the need for decentralization in workflow execution. More recently, several approaches replace a centralized BPEL engine by a set of distributed, loosely coupled and cooperating nodes. It has been identified that one promising solution for this is to use a shared tuplespace for coordination [START_REF] Buhler | Enacting BPEL4WS specified workflows with multiagent systems[END_REF], [START_REF] Martin | A novel approach to decentralized workflow enactment[END_REF], [START_REF] Sonntag | Process space-based scientific workflow enactment[END_REF]. However, in [START_REF] Buhler | Enacting BPEL4WS specified workflows with multiagent systems[END_REF], [START_REF] Martin | A novel approach to decentralized workflow enactment[END_REF], the tuplespace is only used to store data information, our coordination mechanism stores both control and data information into the multiset, which is made possible by the use of the chemical execution model for the coordination of all data and control dependencies. In the scientific area, in [START_REF] Sonntag | Process space-based scientific workflow enactment[END_REF], authors transform a centralized BPEL definition into a set of coordinated processes. Through a shared tuplespace working as a communication infrastructure, the control and data dependencies exchange among processes to make the different nodes interact between them. Again, the use of BPEL hinders from expressing dynamic and selfadaptive behaviors.

As a more general comment, to our knowledge, these work does not provide real software prototype to experimentally validate their proposal.

VII. CONCLUSION

Scientific applications are now built as workflows of services. Workflow management systems gained recently a lot of attention in this context. However, the emergence of new platforms, such as clouds where elasticity and dynamic adaptation are strong requirements, led to a high demand for new models able to represent both wokflows and the platforms, as well as their inherent characteristics.

The chemical model is a promising paradigm naturally capturing parallelism, distribution and dynamics. While the advantages of such a model are now well-established, they still suffer from a lack of proof of concepts and actual deployments.

In this paper, we have proposed a chemistry-inspired workflow management system. A workflow description language and its execution model inspired by such abstractions is discussed. The wide expressiveness (data-flows, control-flows, natural decentralization) of the paradigm is highlighted. Then, its implementation based on the HOCL language, for both centralized and decentralized environment is given. Finally, experiments conducted shown the viability of the concept, lifting a barrier on the path to its actual adoption.

  selectEvens = replace x, ω by ω if x%2 ! = 0 in 1.02 let getM ax = replace x, y by x if x ≥ y selectEvens = s, ω by getM ax, ω 1.07

09 Fig. 1 .

 091 Fig. 1. Chemical workflow representation.

3 .by

 3 01 let invokeServ = replace WSi: CALL:Si, PARAM: in1, . . . , inn , WSi: RESULT:WSi: value , ω 3.04 let preparePass = replace WSi: RESULT:WSi: value , DEST:WSj, ω 3.05 by WSi: PASS:WSj: COMPLETED:WSi: value 3.06 let passInfo = replace WSi: PASS:WSj: ω1 , ω2 , WSj: ω3 3.07 by WSi: ω2 , WSj: ω1, ω3

Fig. 2 .

 2 Fig. 2. Simple merge Pattern.

Algorithm 2

 2 Generic rules -Simple merge pattern 4.01 let sm preparePass = replace DEST:WSj, RESULT:WSi: value 4.02 by PASS:WSj: RESULT:WSi: value , MERGE:Yes 4.03 let sm setFlag = replace-one MERGE:Yes by FLAG INVOKE

Fig. 3 .

 3 Fig. 3. Centralized chemical WMS architecture.

Fig. 4 .

 4 Fig. 4. Decentralized chemical WMS architecture.

Fig. 5 .

 5 Fig. 5. Centralized architecture.

Fig. 6 .

 6 Fig. 6. Decentralized architecture.

Fig. 8 .

 8 Fig. 8. Montage workflow structure

Fig. 7 .

 7 Fig. 7. BlastReport and Cardiac workflows structures

Fig. 9 .

 9 Fig. 9. Performance results, Montage and BlastReport.

Fig. 10 .

 10 Fig. 10. Performance results, CardiacAnalysis.

https://www.irisa.fr/myriads/members/hfernand/hocl/workflows

http://www.myexperiment.org/workflows/2058.html


Samis Trevezas
email: samis.trevezas@ecp.fr

Nikolaos Limnios

Exact MLE and asymptotic properties for nonparametric Semi-Markov models

Keywords: Asymptotic efficiency, Asymptotic normality, Backward recurrence times, Consistency, Maximum likelihood estimation, Nonparametric estimation, Semi-Markov chains AMS Subject Classification: 62F12, 60K15, 62M09

This article concerns maximum likelihood estimation for discrete time homogeneous nonparametric semi-Markov models with finite state space. In particular, we present the exact maximum likelihood estimator of the semi-Markov kernel which governs the evolution of the semi-Markov chain. We study its asymptotic properties in the following cases: (i) for one observed trajectory, when the length of the observation tends to infinity, and (ii) for parallel observations of independent copies of a semi-Markov chain censored at a fixed time, when the number of copies tends to infinity. In both cases, we obtain strong consistency, asymptotic normality, and asymptotic efficiency for every finite dimensional vector of this estimator. Finally, we obtain explicit forms for the covariance matrices of the asymptotic distributions.

Introduction

As opposed to Semi-Markov Processes (SMP) in continuous time, the probabilistic modeling in discrete time has not received analogous attention. One of the earliest references in discrete time is [START_REF] Anselone | Ergodic theory for discrete semi-Markov chains[END_REF], where the author develops some aspects of the ergodic theory of discrete semi-Markov chains (SMCs). When the SMC is not ergodic, but transient with one absorbing state, [START_REF] Gerontidis | Semi-Markov replacement chains[END_REF] studies the long-term behaviour of semi-Markov replacement chains, where independent copies are connected every time an absorption occurs. Applications of semi-Markov models with emphasis on decision processes are given by [START_REF] Howard | Dynamic Probabilistic Systems[END_REF]. Other applications, such as queuing theory and epidemiology, as well as algorithms for implementing SMCs on a computer, are presented by [START_REF] Mode | Computational Methods for Renewal Theory and Semi-Markov Processes with Illustrative Examples[END_REF] and [START_REF] Mode | Stochastic processes in epidemiology[END_REF]. For recent references on theory and applications (e.g. reliability) of discrete SMCs see [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications: Their use in reliability and DNA analysis[END_REF] and [START_REF] Barbu | Discrete time semi-Markov model for reliability and survival analysis[END_REF].

As far as statistical estimation in discrete time is concerned, [START_REF] Barbu | Empirical estimation for discrete time semi-Markov processes with applications in reliability[END_REF] considered empirical estimation for nonparametric semi-Markov models (SMMs) in the case of one observed trajectory. The authors constructed an empirical estimator for the semi-Markov kernel (the kernel that determines the evolution of a SMC) at a finite observation interval, and studied some of its asymptotic properties when the observation length tends to infinity. In particular, they obtained strong consistency and asymptotic normality for each component of the empirical estimator. This estimator can be considered as an approximate maximum likelihood estimator (MLE), since for its construction the sojourn time in the last visited state is systematically neglected. The asymptotic results hold also true for each component of the exact MLE.

In this study, we consider a unified approach in order to tackle the problem of estimation for a single observed trajectory, as well as for multiple trajectories that are generated by independent copies of a semi-Markov chain censored at a fixed time. Both scenarios are interesting for applications. For general references on statistical estimation based on both scenarios, see [START_REF] Kalbfleisch | The statistical analysis of failure time data[END_REF], with emphasis on models and methods for the analysis of failure time data, or [START_REF] Anderson | Statistical models based on counting processes[END_REF]. Let us take an example in reliability. In the case of one observed trajectory, one is mostly concerned with the evolution of an ergodic semi-Markov chain, describing the temporal evolution of one reparable system or one component of the system. The case of L independent copies of a semi-Markov chain arises as well in this context, when one sets into operation L identical components or systems. In this case, there is a useful fixed life time for each system. At the end of this time, the system is rejected independently of its current state. For further applications that concern multiple observed trajectories, see Section 4.

In this article, we present the exact MLE that corresponds to the nonparametric SMM, in order to refine the corresponding estimation of the approximate MLE at a finite interval. We generalise the results of [START_REF] Barbu | Empirical estimation for discrete time semi-Markov processes with applications in reliability[END_REF] by obtaining asymptotic normality for every finite vector of components of the MLE, proving at the same time that the MLE is asymptotically efficient. Additionally, we study the case of parallel observations of independent copies of a SMC and we obtain asymptotic properties, such as strong consistency, asymptotic normality, and asymptotic efficiency, for the part of the semi-Markov kernel that can be estimated until the censoring time. This unified approach reveals that even if the sojourn time in the last visited state can be neglected for obtaining the asymptotic properties of the MLE in the case of one observed trajectory, this is not the case for parallel observations censored at a fixed time. When the last sojourn time is systematically neglected, then the SMCs under study are forced to have a finite support (bounded by the censoring time), and the information for the survival probabilities at the censoring time gets lost.

Our proposed method for recovering this information and obtaining the exact MLE is based on a reparameterisation of the semi-Markov kernel by the parameters that correspond to the natural parametric space of an associated Markov chain (MC), which is formed by the SMC and the sequence of backward recurrence times. The associated MC has an infinite state space under the condition that there exists a state that has conditional sojourn time distribution with an infinite support. Consequently, in this case, the reparameterisation of the semi-Markov kernel by the natural parametric space of the coupled process implies an infinite number of parameters. This justifies the nonparametric framework. The associated MC can play an important role in the understanding of the SMC. Actually, it can be useful, on one hand, to study the probabilistic behavior and limit theorems for SMCs (see, for example, [START_REF] Stenflo | Iterated Function Systems Controlled by a Semi-Markov Chain[END_REF] and Chryssaphinou et al. (2008)), and on the other hand, it can be used to develop statistical inference for the SMMs. The goal of this paper is to develop the latter idea and study MLE for this class of models.

The rest of this paper is organised as follows: In Section 2, we introduce the necessary mathematical notation. In Section 3, we give the exact MLE that corresponds to the semi-Markov kernel of the SMC and we study some of its asymptotic properties. Finally, in Section 4, we study the problem of MLE for parallel observations of independent copies of a SMC censored at a fixed time.

Preliminaries and Notation

Let E be a finite set with cardinality s ∈ N * and M E be the set of real matrices on E ×E. We denote by M E (N) the set of matrix-valued functions defined on N, with values in M E . The following definitions introduce the semi-Markov kernel, the Markov renewal chain (MRC) and the SMC.

Definition 2.1: A matrix valued function q := q(k), k ∈ N ∈ M E (N) is said to be a discrete time semi-Markov kernel if it satisfies the following two properties:

(1)

q ij (k) ≥ 0, i, j ∈ E, k ∈ N, (2) ∞ k=0 j∈E q ij (k) = 1, i ∈ E.
Definition 2.2: Let (Ω, A, P) be a probability space and (J, S) := (J l , S l ) l∈N a bivariate chain with values in E × N. The chain (J, S) is said to be a (homogeneous) Markov renewal chain (MRC) associated to the semi-Markov kernel q if for all k, l ∈ N, i, j ∈ E, it satisfies P-a.s.

P(J l+1 = j, S l+1 -S l = k | J l = i, J l-1 , . . . , J 0 , S l , . . . , S 0 ) = P(J l+1 = j, S l+1 -S l = k | J l = i) = q ij (k).
For the rest we assume that q(0) is the null matrix, and this implies that S is a strictly increasing sequence (therefore excluding instantaneous transitions).

Definition 2.3: The chain Z := (Z n ) n∈N with values in E is said to be a (homogeneous) semi-Markov chain (SMC) associated to (J, S) if

Z n = J N (n) , n ∈ N,
where

N (n) = max{l ∈ N : S l -S 0 ≤ n}, n ∈ N.
(1)

The process S keeps track of the successive time points that changes of states in Z occur (jump times), and J records the visited states at these time points. We also make the assumption that S 0 = 0, and this means that the first observation of the SMC is considered to be on a jump time. The sequence defined by (1) is the counting process of the number of jumps (changes of states) until time n. Let also X := (X l ) l∈N * , where X l+1 = S l+1 -S l . The process X corresponds to the sequence of total sojourn times in the successively visited states.

Remark 1 : We want a jump time to indicate a change of state, so we will assume that q ii (k) = 0, for all i ∈ E, k ∈ N. This assumption excludes absorbing states from our analysis, but in any case, we can modify easily our setting by partitioning the state space in a subset where this condition holds true and in the complementary set, where only absorbing states are included and satisfy k q ii (k) = 1. In this case, each absorbing state forms a simple renewal chain.

For the statistical analysis, we are interested in the nonparametric semi-Markov model, where the SMC Z is defined on the probability space (Ω, A , P q ), and the semi-Markov kernel q is identified with the infinite dimensional vector q := q ij (k) i,j∈E,k∈N . We also assume that its components keep the lexicographical order of the triplet (i, k, j).

Let E -i := E \ {i}. We denote by

Q i :=    q i = q ij (k) ∈ R (N * ) s-1 ; q ij (k) ≥ 0, j∈E-i k∈N * q ij (k) = 1    ,
for all i ∈ E, the parametric space of all the possible discrete distributions on the set (N * ) s-1 , where the components of the vector q i are considered to be ordered with respect to the lexicographical order of (k, j) ∈ N * × E. Then, by taking s independent copies of the space Q i , we will represent the semi-Markov kernel q as a point in the parametric space Q := i∈E Q i .

Let U := (U n ) n∈N be the sequence of backward recurrence times (BRT), where

U n = n -S N (n) .
Note that U 0 = 0, since S 0 = 0 by assumption. Each BRT U n keeps track of the elapsed time after the last jump of the process Z until time n. The pair (Z, U) can be shown to be a Markov chain (e.g., [START_REF] Limnios | Semi-Markov Processes and Reliability[END_REF], Theorem 3.12). See also Chryssaphinou et al. (2008), where the authors study some properties of (Z, U). In this paper, the identically non-zero transition probabilities of the MC (Z, U) will serve as a suitable reparameterisation of the semi-Markov kernel q, in order to obtain its exact MLE. Now, we denote by Hi (•) the survival function in state i, which is given as a function of the semi-Markov kernel as follows:

Hi (u) := P q (X l+1 > u | J l = i) = 1 - j∈E u k=0 q ij (k), u, l ∈ N.
(2)

If (Z, U) is irreducible and positive recurrent, then the corresponding stationary proba-bilities are given by

π i,u = Hi (u) µ ii , (i, u) ∈ E × N,
where µ ii corresponds to the mean recurrence time between two successive visits of the MC (Z, U) to the state (i, 0). The mean recurrence time has the following representation

µ ii = E(G i | Z 0 = i, U 0 = 0), (3)
where G i is the first hitting time of state i ∈ E, after leaving the initial state, that is,

G i = inf{n ≥ S 1 : Z n = i}.
For the remainder let T m := {0, 1, . . . , m-1} and T + m be the set T m translated one unit to the right. Let also E 2 * := E 2 \ {(i, i), i ∈ E}. Now, we define some counting processes that will be used in the sequel, appropriate for the pair (J, X) :

N ij (k, m) := N (m)-1 l=0 1 {Jl=i,Jl+1=j,Xl+1=k} , (i, j) ∈ E 2 * , k ∈ T + m , (4)
N i (m) := N (m)-1 l=0 1 {Jl=i} . i ∈ E, (5)
The process N ij (k, m) counts the number of (one-step) transitions of J, from i to j, that took place in k time units, until time m. The process N i (m) counts the number of visits of J to state i, until time m, not counting the last visit J N (m) . Note that the time step is given by the process Z.

Similarly, we define some counting processes, appropriate for the pair (Z, U) :

N B i,u (j; m) := m-1 n=0 1 {Zn=i,Un=u,Zn+1=j} , (i, j) ∈ E 2 , u ∈ T m , (6)
N B i,u (m) := m-1 n=0 1 {Zn=i,Un=u} , (i, u) ∈ E × T m . (7)
The process N B i,u (j; m) counts the number of (one-step) transitions of (Z, U), from (i, u) to (j, k), until time m. Since k is uniquely determined when the next j is known, we refer only to j for notational convenience. The process N B i,u (m) counts the number of visits of (Z, U) to state (i, u), until time m -1. We note that the superindex B (from Backward) is only used in order to distinguish easier the counting processes that refer to (Z, U).

If we denote by

M (s) := inf{m ∈ N * : N i (m) > 0, ∀i ∈ E},
then for m ≥ M (s) , every state and at least one transition from every state of the MC J to some other state have been recorded.

Exact MLE based on a single trajectory

In this Section, we construct the exact MLE of the semi-markov kernel, when a fixed length observation Z m 0 := (Z 0 , Z 1 , . . . , Z m ) of a SMC Z is available, and we study its asymptotic properties as m tends to infinity.

The likelihood function L m (q) := P q (Z m 0 ) that corresponds to an observation Z m 0 can be written as:

L m (q) = a(Z 0 )   N (m) k=1 q Jk-1Jk (X k )   HJN(m) (U m ), (8)
where a(•) is the initial distribution and Hi (•) defined as in Relation (2).

Since the factor HJN(m) (U m ) that appears in equation ( 8) corresponds to the sojourn time in the last visited state of the system until time m, the information that we get for the semi-Markov kernel from this part can be neglected for large m for positive recurrent SMCs (see [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications: Their use in reliability and DNA analysis[END_REF]). This leads to an approximate MLE (quasi-MLE), which is asymptotically equivalent with the exact one in the case of a single trajectory. [START_REF] Barbu | Empirical estimation for discrete time semi-Markov processes with applications in reliability[END_REF] used this fact and considered the approximate MLE in order to develop statistical inference for nonparametric SMMs. Under the assumption that the initial distribution does not depend on q, and by ( 4), the authors reduced the maximisation problem to the maximisation of an approximate log-likelihood function:

lm (q) = i,j∈E 2 * k∈T + m N ij (k, m) log q ij (k). (9)
In this way, they obtained the following M-estimator (see [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF] for a general treatment of M-estimators):

qij (k, m) = N ij (k, m) N i (m) . ( 10
)
If N i (m) = 0, we take by convention 0/0 = 0.

Remark 1 : This convention allows us to extend the estimation to the case that we have not observed all the states (m < M (s) ). This is valid if we assume that the considered s-state SMCs that we deal with, are inside a larger parametric space that includes all the SMCs with cardinality bounded by s, or equivalently, the cardinality of the state space of the considered SMCs is a parameter bounded by s. This remark extends and proves to be fruitful even in the case of the chain (Z, U). This pair forms an infinite state space Markov chain on condition that there exists at least one state of the SMC Z with unbounded sojourn time. In this case, the state space of the considered MCs should be considered as a nonnegative integer parameter (not bounded).

Construction of the MLE

In this subsection, we infer the exact MLE of this model and its relationship with the M-estimator (empirical estimator or approximate MLE) given by ( 9).

Let

p i,u;j := P q (Z n+1 = j, U n+1 = 0 | Z n = i, U n = u) if j = i, P q (Z n+1 = i, U n+1 = u + 1 | Z n = i, U n = u) if j = i.
(11)

We will reparameterise the model by the natural parameter that corresponds to the MC (Z, U), that is, η := (p i,u;j ) i,j∈E,u∈N , where all the identically zero probabilities of the transition matrix of (Z, U) have been excluded. Note also that q ≡ Φ(η), where the component functions of Φ for i = j are given by

Φ ijk (η) = p i,0;j if k = 1, p i,k-1;j k-2 u=0 p i,u;i if k ≥ 2. ( 12
)
Each one of the functions Φ ijk depends actually on k variables (finite number). Let η m := (p i,u;j ) i,j∈E,u∈Tm , where η m parameterise the first sm distributions with finite support (see Relation ( 11)), for which we can have information until time m.

In the following proposition, we will obtain the MLE for the above distributions and we will transfer the results to the MLE for the parameters q m := q ij (k) (i,j)∈E 2 * ,k∈T + m of the semi-Markov kernel q that can be estimated until time m. Additionally, we connect the empirical estimators qij (k, m) given by ( 10) with the MLE q ij (k, m) obtained here. Proposition 3.1: i) The MLE that corresponds to a fixed length observation Z m 0 of a nonparametric SMM, driven by the semi-Markov kernel q, is given for i = j by

q ij (k, m) = Φ ijk ( η m (m)) = p i,k-1;j (m) H i (k -1, m), k ∈ T + m , (13)
where

H i (k -1, m) = k-2 u=0 p i,u;i (m).
ii) The exact relationship between the MLE q ij (k, m), given by ( 13), and the approximate MLE qij (k, m), given by (10), reads for (i, j, k) such that qij (k, m) > 0 as follows

q ij (k, m) = qij (k, m) 1 - 1 {Zm=i,Um>0} N i (m) + 1 1 + 1 {Zm=i,k>Um>0} N B i,Um (m) . ( 14
)
Proof : By using the counting processes (6), the likelihood function L m (q) given by ( 8) can be written as a function of η m in the form

L m (η m ) = a(Z 0 ) i,j∈E u∈Tm p N B i,u (j;m) i,u;j .
Since the initial distribution does not vary with η m , the log-likelihood function that has to be maximized is given by

l m (η m ) = i,j∈E u∈Tm N B i,u (j; m) log p i,u;j .
(15) S. Trevezas and N. Limnios

The dependence relations that we have to take into consideration in a nonparametric framework are imposed from the stochastic nature of the transition probability matrix P of the MC (Z, U). For each i ∈ E, u ∈ N, j∈E p i,u;j = 1. ( 16)

By using the expressions ( 6) and ( 7), then by maximizing with respect to η m , under the constraints ( 16), the likelihood function l m (η m ) given by ( 15), we get easily

η m (m) = p i,u;j (m) = N B i,u (j; m) N B i,u (m)
.

From the invariance principle of MLE and Relation (12) we conclude (13). ii) By using ( 6) and ( 7), Relation ( 13) can be written in the form

q ij (k, m) = N B i,k-1 (j; m) N B i,k-1 (m) N B i,k-2 (i; m) N B i,k-2 (m) • • • N B i,1 (i; m) N B i,1 (m) N B i,0 (i; m) N B i,0 (m) , k ≥ 2. ( 17
)
Notice that

N B i,k-1 (j; m) = N ij (k, m), (18)
since in order to have a transition from i to j in k time units for the process Z, a transition from (i, k -1) to (j, 0) has to take place for the process (Z, U). Moreover, for 1 ≤ l ≤ k -1 we have

N B i,l-1 (i; m) = N B i,l (m) + 1 {Zm=i,Um=l} , (19)
and finally

N B i,0 (m) = N i (m) + 1 {Zm=i,Um>0} . (20)
By inserting ( 18)-(20) in Relation (17) we get

q ij (k, m) = N ij (k, m) N B i,Um (m) + 1 {Zm=i,k>Um>0} N i (m) + 1 {Zm=i,Um>0} N B i,Um (m) . ( 21
)
By equation ( 21) we conclude easily (14).

Remark 2 : Note that N B i,Um (m) = 0 iff at time m the system stayed for the first time at state i for U m time units. Consequently, the subsequence of m that satisfies this equation defines the time points that a new state of (Z, U) has been recorded. Also, at these time points we have that 1 {Zm=i,k>Um>0} = 0 for k such that qij (k, m) > 0, since at time m no transition has occurred from state i to state j in more than U m time units. Therefore, the right-hand member of equation ( 14) is well defined by the convention 0/0 = 0.

Asymptotic properties of the MLE

In this subsection, we study some asymptotic properties of the nonparametric MLE of the semi-Markov kernel. Theorem 3.4 summarises the asymptotic results in the framework of convergence of stochastic processes, since we have an infinite dimensional parametric space. Under a suitable condition, the MLE is strongly consistent and every finitedimensional subvector of the MLE is asymptotically normally distributed and asymptotically efficient. In order to prove asymptotic normality, we use the reparameterisation (12) of the semi-Markov kernel. For every fixed ñ ∈ N, we obtain in Lemma 3.2 the asymptotic distribution of the MLE of the vector parameter η ñ and then we transfer this result to the MLE of the vector parameter q ñ (Proposition 3.3). Now, we state the following condition concerning the subclass of SMCs to be considered: Condition (A1): The SMC Z under P q is irreducible and positive recurrent.

The above condition is necessary in the case that we search the long term behavior of the MLE for one trajectory of the SMC Z. This condition will be relaxed in the next Section, where we study the asymptotics of the MLE for an increasing number of independent realisations of a SMC Z, censored at a fixed time. Also, for notational simplicity, we assume in this Section that all the elements of the semi-Markov kernel are strictly positive under the true value. This assumption implies that condition (A1) holds. Otherwise, the following Lemma holds true for the MLE that corresponds to the probabilities that are strictly positive under the true value. The proof of this Lemma is deferred to Appendix A.

Lemma 3.2: Let (Z, U) be the MC associated with the SMC Z, where U is the sequence of BRT. Under condition (A1), the sequence of random vectors ξ ñ(m) := √ m( η ñ(m)η ñ), converges P η -weakly, as m tends to infinity, to a random vector ξ ñ, where ξ ñ ∼ N (0 , Γ ñ) and Γ ñ is a covariance matrix of dimension (s 2 ñ) × (s 2 ñ), defined in a block diagonal form

Γ ñ = diag{ 1 π i,u Λ iu ; (i, u) ∈ E × T ñ}, (22)
where the matrix Λ iu , for each (i, u) ∈ E × T ñ is given by

Λ iu = p i,u;j [δ jl -p i,u;l ] j,l∈E . (23)
The proof of the following Proposition is based on Lemma 3.2 and it is deferred to Appendix A.

Proposition 3.3: Under condition (A1), the sequence of random vectors ξ ñ(m) := √ m( q ñ(m)q ñ), converges P q -weakly, as m tends to infinity, to a random vector ξ ñ, where ξ ñ ∼ N (0 , ∆ ñ) and ∆ ñ is a covariance matrix of dimension (s 2s)ñ × (s 2s)ñ, defined in a block diagonal form

∆ ñ = diag µ ii K i ; i ∈ E , ( 24
)
where the matrix K i , for each i ∈ E is given by

K i = q ij (k)[δ (j,k)(r,l) -q ir (l)] (j,k),(r,l)∈E-i×T + ñ . ( 25
)
and µ ii is given by (3).

Remark 3 : A direct proof of Proposition 3.3 is also possible by means of the technique that we use in order to prove Lemma 3.2. In this case, the martingale CLT can be applied by considering an appropriate martingale difference with respect to the natural filtration generated by the Markov chain (J, X). A possible choice of the martingale difference (compare with (A8), in the case of the MC (Z, U)) is given as follows:

Y n := 1 {Xn≤ñ} [t Jn-1JnXn -E q (1 {Xn≤ñ} t Jn-1JnXn | J n-1 )]. ( 26
)
Nevertheless, our method gives a unified treatment of the study of the MLE, when we have a single or multiple trajectories.

In order to prove asymptotic efficiency, we will use the same reparameterisation (12). Since this property is parameterisation independent, it is enough to show that for an arbitrary ñ ∈ N the MLE of the vector parameter η ñ is asymptotically efficient. This will be true, if the asymptotic covariance matrix of the MLE of the vector parameter η ñ is the inverse of the asymptotic FIM for this model. Then, the MLE of the vector parameter q ñ will be asymptotically efficient as well.

Let η * m be a minimal representation of η m , by excluding for each (i, u) the parameter p i,u;i (see Relation ( 16), for the dependence relations). In order to compute the asymptotic Fisher Information Matrix (FIM), we define first the following matrix:

I(η * m ) := E α,η -∂ 2 ℓ m (η m ) ∂p i1,u1;j1 ∂p i2,u2;j2 (i1,u1,j1),(i2,u2,j2) , ( 27
)
where the admissible indices (i 1 , u 1 , j 1 ) and (i 2 , u 2 , j 2 ) are these that correspond to η * m , and ℓ m (η m ) is given by (15). The matrix I(η * m ) is an unnormalised FIM that corresponds to the trajectory Z m 0 , when the scenario of a single trajectory of a SMC is considered. When multiple trajectories (i.i.d.) of a SMC censored at time m are considered (see Section 4), then it is the exact FIM that corresponds to this scenario. The matrix

I m (η * ñ) := m -1 I(η * ñ) (28)
is a normalised version (information per unit time) of the matrix defined by ( 27), where now the admissible indices (i 1 , u 1 , j 1 ) and (i 2 , u 2 , j 2 ) are these that correspond to η * ñ (recall that η * ñ parameterises the first sñ transition probability distributions of (Z, U)). Now, let for m ∈ N * , q(m) = q m (m), 0, 0, . . . , that is, the sequence of estimators that can be decomposed into two parts, a finite part that corresponds to the MLE of the vector q m , and the second part that completes with infinite zeros the parameters that we have no information up to time m. Theorem 3.4 : Under condition (A1), i) the sequence of MLE q(m) is a strongly consistent estimator of q, that is, as m → ∞, q(m) Pq-a.s.

----→ q, ii) the sequence of random processes ξ(m) := ξ m (m), 0, 0, . . . , m ∈ N * , where ξ m (m) := √ m( q m (m)q m ), converges P q -weakly, as m tends to infinity, to a random process ξ, where for any ñ ∈ N * , the corresponding random vector ξ ñ ∼ N (0 , ∆ ñ), and ∆ ñ is given by ( 24), iii) every finite-dimensional subvector of the MLE q(m) is asymptotically efficient.

Proof : i) It is straightforward to see by Relation ( 14) that the MLE q ij (k, m) and the empirical estimator qij (k, m) of the parameter q ij (k, m) converge P q -a.s. to the same limit, if that limit exists. [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications: Their use in reliability and DNA analysis[END_REF] have proved that qij (k, m) is strongly consistent (see Theorem 4.1.). Therefore, the strong consistency of q(m) follows from the coordinatewise consistency.

ii) The sequence of random processes ξ(m) converges P q -weakly to the process ξ iff for any ñ ∈ N * the sequence of random vectors ξ ñ(m) converges P q -weakly, as m tends to infinity, to the corresponding finite dimensional vector ξ ñ. This result is a consequence of Proposition 3.3. iii) Since the asymptotic efficiency of the MLE is parameterization independent, it suffices to show that the MLE η ñ(m) of the vector parameter η ñ (for an arbitrary ñ) is asymptotically efficient. This is true if the limiting matrix

I ∞ (η * ñ) := lim m→∞ I m (η * ñ),
where I m (η * ñ) is given by ( 28), exists and its inverse equals the asymptotic covariance matrix Γ ñ given by ( 22). The proof of this assertion is deferred to Appendix A.

Independent observations of a SMC censored at a fixed time

In this Section, we consider L trajectories of independent copies of a SMC Z censored at a fixed time m. The setting of multiple independent trajectories arises in many situations, and statistical inference should be developed in this framework, when each trajectory is assumed to follow semi-Markov dynamics. We give here some relative examples: i) in reliability, when L identical but independent systems are observed for a fixed time period (see Introduction), ii) in epidemiology (biomedical interest), when L patients, victims of a certain disease, are observed during disease progression (see, e.g., [START_REF] Commenges | Multi-state Models in Epidemiology[END_REF], for the so-called multistate models). Nevertheless, in this case, the SMCs are absorbing, including a fatal state, iii) In DNA analysis, when L large sets of trusted examples (called training sets and correspond here to fixed-length trajectories) are typically used in order to estimate the parameters of a probabilistic model that underlies the generation of a sequence of DNA bases. Hidden Markov models and hidden semi-Markov models have been used extensively in this context (see [START_REF] Durbin | Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids[END_REF] and [START_REF] Burge | On Discrete Time Semi-Markov Chains and Applications in Words occurrences[END_REF]). In the case of hidden semi-Markov models, there is a hidden semi-Markov chain that models the succession of different regimes (evolutionary states), where the emission probability distribution of the DNA bases changes. MLE based methods can be used in order to combine the different training sets and obtain the parameter estimates, iv) in Monte-Carlo EM (Expectation-Maximisation) type algorithms for hidden semi-Markov models (intended for inference with a large data set), where independent trajectories of a SMC are generated (simulated) at each iteration of the EM algorithm (at each E-step, and under the current value of the parameter), in order to approximate the so-called Qfunction of the EM-algorithm (equivalent to the incomplete likelihood function, as far as maximisation is concerned) that corresponds to this model (see, e.g., [START_REF] Malefaki | An EM and a Stochastic Version of the EM Algorithm for Nonparametric Hidden Semi-Markov Models[END_REF]).

In the sequel, we study the MLE that corresponds to this model and we are interested in its asymptotic properties. We denote by (Z m 0 ) l := (Z

(l) 0 , Z (l) 1 , . . . , Z (l)
m ) the l-th realisation of independent copies of a SMC Z censored at the fixed time m. If the problem was to make inferences for a class of non positive recurrent SMCs, then the estimation based on a single trajectory of a SMC is not satisfactory. On the other hand, when many trajectories corresponding to independent copies of a CSM are available, the structure of the model can be better explored. For this reason, we make no assumption here for the type of communication among the different states of the SMC Z. For state classification of a finite-state SMC and different types of communication, see, e.g., Definition 3.11 in [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications: Their use in reliability and DNA analysis[END_REF]. Now, we give a suitable condition for this type of estimation.

Condition (A2):

The initial distribution α(•) under P α,q i) has a support that gives strictly positive probability to every communication class of Z, ii) does not depend on q.

Interpretation. Since in each trajectory, it is only the first observed state that determines the set of observed states, it suffices to assume that at least one state for each communication class of Z is accessible through the initial observation. This is done in part i) of the above assumption, which guarantees that all classes will be observed (with probability one) after a sufficiently large L . Besides, due to part ii) we can decompose the maximisation problem into two independent parts. The first part corresponds to the initial distribution and the other part to the semi-Markov kernel.

Despite the advantage of the arbitrary structure of the system, nothing guarantees that the complete structure of the system will be revealed even after a great number of observed trajectories. The censoring time m plays an important role on the transitions that can actually be observed. Naturally, we cannot observe transitions from a state i to a state j when the minimum time of doing so, exceeds m, even if j is accessible from i under the true value of the parameter. Consequently, our target is to estimate the part of the semi-Markov kernel that is accessible until time m. By assuming that S 0 = 0, this part corresponds to q m . We prove in the sequel (see Theorem 4.4) that the MLE of the parameter q m is strongly consistent, asymptotically normal and asymptotically efficient.

Additionally, this type of estimation allows to measure the concentration of the sojourn time distribution for each i ∈ E on the finite set {1, 2, . . . , m}. This can be achieved by the MLE of the distribution function of the sojourn time in each state at time m. The subsequent asymptotic normality of the MLE (see Corollary 4.5) enables us to construct asymptotic confidence intervals in order to measure the error of this estimation.

In this Section, we indicate how we can treat generally the case where under the true value of the semi-Markov kernel not necessarily all the elements are strictly positive. To this end, we denote by

D * i|m := {(j, k) ∈ E × T + m : q ij (k) > 0} and d * i|m := |D * i|m |.
Moreover, the part of q m that corresponds to the strictly positive parameters q ij (k) is denoted by q + m . With the same reasoning the vector η + m denotes the probabilities of the MC (Z, U) that are strictly positive. Additionally, η * m designates η + m in a minimal representation by excluding for each (i, u) the parameter p i,u;i .

Construction of the MLE

By following the lines of the previous Section we use the reparameterisation of the SMM by the natural parametric space of the MC (Z, U) censored at time m. Let N B,l i,u and N B,l i,u;j be the l-th realisations of the counting processes N B i,u and N B i,u;j respectively, defined in the previous Section by ( 7) and ( 6). Moreover, we define the counting processes

N B i,u (L) := L l=1 N B,l i,u , where N B i,u (m; L) = L l=1 N B,l i,u (m), (29)
and

N B i,u;j (L) := L l=1 N B,l i,u;j , where N B i,u (j, m; L) = L l=1 N B,l i,u (j, m). (30)
The r.v. N B i,u (m; L) indicates the number of visits at state (i, u) ∈ E × N, until time m ∈ N * for the L trajectories, by neglecting the last visited state (Z m , U m ). Besides, the r.v. N B i,u (j, m; L) indicates the number of possible transitions of the type (i, u), (j, k) (that is, (j, k) = (j, 0), j = i, or (j, k) = (i, u + 1)), until time m for the L trajectories.

This type of estimation enables us to estimate the initial distribution α(•) of the SMC Z. To this end, we define for each i ∈ E the counting process N α i := {N α i (L)} L∈N * of the times that i appeared as initial state in the L trajectories, where

N α i (L) = L l=1 1 {Z (l) 0 =i} . (31)
The following proposition extends the results of Proposition 3.1 for MLE from 1 to L trajectories of a SMC Z censored at time m. Moreover, it is not necessary to assume that Z is positive recurrent since by the independence of the copies of Z the initial states can belong to different communication classes.

Proposition 4.1: Let Z be a SMC defined on (Ω, F, P α,q ) (nonparametric model), with state space E = {1, 2, . . . , s}. Under condition (A2), the MLE ( α(L), q m (m; L)) that corresponds to L trajectories (Z m 0 ) i for i ∈ T + L from L independent copies of a SMC Z censored at time m, is given by its components, for (i, j) ∈ E 2 * , k ∈ T + m , as follows

q ij (k, m; L) = Φ ijk ( η m (m; L)) = p i,k-1;j (m; L) H i (k -1; L), (32)
α i (L) = N α i (L) L , (33)
where

H i (k -1; L) = k-2 u=0 p i,u;i (m; L)
and the MLE η m (m; L) of the parameter η m is given by

η m (m; L) = p i,u;j (m; L) = N B i,u (j, m; L) N B i,u (m; L) . ( 34
)
If N B i,u (m; L) = 0, then by convention 0/0 = 0.

Proof : The likelihood function L L (α, q) for L independent observations (Z m 0 ) i , i ∈ T + L , is given by L L (α, q) = L l=1 α J (l) 0   N (l) (m) k=1 q J (l) k-1 J (l) k (X (l) k )   HJ N (l) (m) (U (l) m ).
By ( 31) and ( 30), the likelihood function L L (α, q) is given as a function of (α, η m ) by

L L (α, η m ) = i∈E α N α i (L) i i,j∈E l∈T + L u∈Tm p N B,l i,u (j,m) i,u;j .
Consequently, we obtain the following log-likelihood function

ℓ L (α, η m ) = i∈E N α i (L) log α i + i,j∈E l∈T + m u∈Tm N B,l i,u (j, m) log p i,u;j . (35)
The dependence relations that we have to take into consideration in a nonparametric framework are given by ( 16) for η m and similarly for the initial vector α. From these dependence relations, condition (A2), ( 29) and by maximising the log-likelihood function (35) with respect to α and η m , we obtain (33) and (34). From the invariance principle of the MLE and Relation (12) we infer (32).

Asymptotic properties of the MLE

In this subsection, we study some asymptotic properties of the MLE of the parameter (α, q m ). In particular, we prove strong consistency, asymptotic normality and asymptotic efficiency of the MLE. These properties will be deduced from the corresponding properties of the MLE of the parameter (α, η m ), which we present in the following lemma.

Remark 1 : Under condition (A2), the MLE α(L) is strongly consistent, asymptotically normal and asymptotically efficient. Strong consistency is a direct consequence of the strong law of large numbers, applied to the sequence of i.i.d. random vectors {(1 {Z (l) 0 =i} : i ∈ E)} l∈N * . Asymptotic normality follows from the CLT for i.i.d. random vectors, that is, the vector ζ(L) := √ L( α(L)α) converges P α,η -weakly, when L → ∞, to ζ ∼ N (0 , Γ (L) ), where Γ

(L) ij := Cov α,η (ζ i , ζ j ) = α i (δ ij -α j )
. Asymptotic efficiency also follows from standard arguments for the MLE, in the i.i.d. framework.

From the above remark, and since the MLE α(L) and the MLE η + m (m; L) are asymptotically independent, we will ignore in the rest the parameter α from the estimation, and we will focus our attention on the parameter η m (or η + m ). The part ii) of the lemma is an adaptation of the theorem of [START_REF] Anderson | Statistical inference about Markov Chains[END_REF] to our case. For the rest, we denote by d i,u the number of different j such that p i,u;j is strictly positive. m is a covariance matrix given in a block diagonal form as follows

Γ (L) m = diag{ 1 π i,u|m Λ iu : (i, u) ∈ E × T m }, (36)
where for each (i, u) ∈ E × T m , the matrix Λ iu of dimension d i,u × d i,u is given by ( 23) and π i,u|m is given by

π i,u|m := E α,η N B i,u (m) = r∈E n∈Tm α r p (n) r,0;i,u , (37)
where

p (n) r,0;i,u = P α,η (Z n , U n ) = (i, u)|(Z 0 , U 0 ) = (r, 0)),
iii) the MLE η * m (m; L) of the parameter η * m is asymptotically efficient. Proof : i) Notice by the form of η m (m; L), that if we divide each term of the fraction with L, and apply the strong law of large numbers, we have as L → ∞,

p i,u;j (m; L) Pα,η-a.s. -----→ E α,η N B,l i,u (j, m) E α,η N B,l i,u (m) , (38)
since both counting processes are sequences of i.i.d. bounded (by m) r.v. . By (6), we have also that

E α,η N B,l i,u (j, m) = p i,u;j m-1 n=0 P α,η (Z (l) n , U (l) n ) = (i, u) = p i,u;j E α,η N B,l i,u (m). ( 39
)
By inserting (39) in (38), we obtain the result. ii) [START_REF] Anderson | Statistical inference about Markov Chains[END_REF] proved the asymptotic normality of the MLE of a finite MC even in the inhomogeneous case. At first, they considered known the vector of observed frequencies for the initial states in L trajectories. Next, they proved the S. Trevezas and N. Limnios asymptotic results when this vector is multinomially distributed. The steps of their proof can be applied to the case of L trajectories corresponding to L independent copies of the MC (Z, U) observed until time m, since they are based on the multidimensional CLT for i.i.d. random vectors and not on the structure of (Z, U). Nevertheless, in order to apply them in our case, it only suffices to note that the support of the initial distribution reduces to the states {(i, 0)} i∈E .

iii) Since we are in an iid setting, in order to compute the FIM, it is enough to take the log-likelihood contribution of a single term, corresponding here to the vector Z m 0 . Therefore, the FIM is the matrix I(η * m ), defined by ( 27). We have directly that

∂ℓ m (η m ) ∂p i1,u1;j1 = N B i1,u1 (j 1 , m) p i1,u1;j1 - N B i1,u1 (i 1 , m) p i1,u1;i1 ,
and also that

∂ 2 ℓ m (η m ) ∂p i1,u1;j1 ∂p i2,u2;j2 =              - N B i1,u1 (j 1 , m) p 2 i1,u1;j1 - N B i1,u1 (i 1 , m) p 2 i1,u1;i1 , if (i 1 , u 1 , j 1 ) = (i 2 , u 2 , j 2 ), - N B i1,u1 (i 1 , m) p 2 i1,u1;i1 , if (i 1 , u 1 ) = (i 2 , u 2 ), j 1 = j 2 , 0, if (i 1 , u 1 ) = (i 2 , u 2 ).
(40) By taking the negative expectations of the terms (40) and by ( 39) and ( 37), we deduce that

I(η * m ) = diag{I i,u (η * m ); (i, u) ∈ E × T m }, (41)
where

I i,u (η * m ) := π i,u|m p -1 i,u;i + δ j1j2 p -1 i,u;j1 j1,j2 . (42)
By (36), the asymptotic covariance matrix Γ (L) m of the MLE, considered in a minimal representation (by excluding the lines and columns that correspond to p i,u;i ), is the inverse of the FIM I(η * m ) given by ( 41) and ( 42), iff the matrix Λ iu given by ( 23), is the inverse of π -1 i,u|m I i,u (η * m ) given by ( 42). This can be verified with elementary calculus. The following definition introduces the Markov renewal function that plays a very important role in the theory of Markov renewal chains. This definition was given by [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications: Their use in reliability and DNA analysis[END_REF] in a different context, in the study of one trajectory.

Definition 4.3: The matrix function Ψ = Ψ(m), m ∈ N ∈ M E (N), where Ψ ij (m) := E(N B j,0 (m + 1) | Z 0 = i)
is called the discrete time Markov renewal function. We also denote by

Ψ i (m) := EN B i,0 (m + 1). ( 43
)
With the following theorem we adapt the asymptotic results of the MLE for the MCs to the case of SMCs by making also explicit the covariance matrix form of the asymptotically normal distribution of the MLE. The trick consists in showing that this theorem can be proved by Lemma 4.2 in a similar way that Proposition 3.3 can be proved by Lemma 3.2.

Theorem 4.4 : Under condition (A2), i) the MLE q m (m; L) of the parameter q m , given by (32), is strongly consistent, ii) the sequence of random vectors ξ m (L) := √ L( q + m (m; L)q + m ), converges P α,η -weekly, when L → ∞, to a random vector ξ m , where ξ m ∼ N (0 , ∆

L) m ) and ∆ (L) (
m is a covariance matrix given in a block diagonal form as follows

∆ (L) m = diag 1 Ψ i (m -1) K i : i ∈ E , (44)
where for each i ∈ E the matrix K i of dimension d * i|m ×d * i|m is given by ( 25) and Ψ i (m-1) is given by ( 43) and assumed to be strictly positive. iii) the MLE q + m (m; L) of the parameter q + m , is asymptotically efficient. Proof : i) By Lemma 4.2, we have that the MLE η m (m; L) of the parameter η m , given by (34), is strongly consistent. Let Φ m be the restriction of Φ, given by ( 12), to the domain that corresponds to the vector parameter η m in a minimal representation. By applying the continuous mapping theorem for the function Φ m we obtain the desired result. ii) Notice by ( 37) and ( 43) that

π i,0|m = Ψ i (m -1). (45)
Moreover,

E α,η N B i,u (m) = m-1 n=0 P α,η (Z n = i, U n = u) = m-1 n=u P α,η (Z n-u = i, U n-u = 0) × P α,η (X N (n-u)+1 > u|Z n-u = i, U n-u = 0) = Hi (u) m-1 n=0 P α,η (Z n = i, U n = 0) = Hi (u) E α,η N B i,0 (m).
This implies by ( 37) and ( 45) that

π i,u|m = Hi (u)Ψ i (m -1). ( 46
)
We deduce that the ratio of Hi (u) with π i,u|m does not depend on u. By inserting

(m, π i,u|m , (Ψ i (m -1)) -1 , L, Γ (L)
m ) in the place of (ñ, π i,u , µ ii , m, Γ ñ) in the statement of Lemma 3.2, the steps of the proof needed for the part ii) of this theorem can be deduced by Lemma 4.2, in the same way that Proposition 3.3 can be deduced by Lemma 3.2. iii) The proof is a direct consequence of the part iii) of Lemma 4.2, since the asymptotic efficiency of the MLE is parameterisation independent.

Remark 2 : i) The term Ψ i (m -1) that appears in the expression (44) and which is defined by (43), can be expressed as a function of the semi-Markov kernel q. If q (n) represents the n-th convolution power of q, (see [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications: Their use in reliability and DNA analysis[END_REF], pages 48-49), then the following relation holds

Ψ i (m -1) = r∈E n,l∈Tm α r q (n) ri (l). (47)
Consequently, by ( 47) and ( 46) we obtain an alternative form of (37),

π i,u|m = Hi (u) r∈E n,l∈Tm α r q (n) ri (l).
ii) The term Ψ i (m-1) by ( 43) is strictly positive iff there exists a state j where α j > 0 and the state i ∈ E is accessible from j until time m. If the support of the initial distribution is E, this is always true. Otherwise, in Theorem 4.4 we have to consider only the states that satisfy the above property.

In order to measure the concentration of the sojourn time distribution at each state i ∈ E inside the set {1, 2, . . . , m}, we consider the MLE of the distribution function H i (m) at time m. Since

H i (m) = (j,k)∈D * i|m q ij (k), ( 48
)
the MLE is given by H i (m; L) = (j,k)∈D * i|m q ij (k, m; L).

This additive form is convenient in order to deduce easily the asymptotic results of strong consistency and asymptotic normality of the MLE. The proof is deferred to Appendix A.

Z (L) m = diag H i (m) Hi (m) Ψ i (m -1) : i ∈ E . ( 49
)
allows us to apply the Cramér-Wold device, where we need to show that for an arbitrary t ∈ R (s 2 -s)ñ , the r.v. m -1/2 tζ * (m) converges in distribution to the r.v. tζ * . For this purpose, we will represent the sequence of r.v. {tζ * (m)} m≥1 as the sequence of partial sums of another sequence {Y n } n≥1 that has the property to be a martingale difference adapted to the natural filtration of the underlying probability space (Ω, A , P η ) of the Markov chain (Z, U). Let F n := σ{(Z k , U k ) 0≤k≤n }, that is, the σ-algebra generated by the random vector (Z k , U k ) 0≤k≤n . For each n ∈ N * we define the random variable Let also t iu := t iuj ; j ∈ E -i , and Λ * iu ∈ R (s-1)×(s-1) the matrix that results from Λ iu , which appears in Relation (22) if we extract the i-row and the i-column. Then, the coefficients λ i,u (t, η) given by (A15) can be expressed as follows λ i,u (t, η) = t iu Λ * iu t ⊤ iu .

It is easy to see that det(Λ * iu ) = j∈E-i p i,u;j 1 -j∈E-i p i,u;j = j∈E p i,u;j .

Therefore, since Λ * iu is a symmetric matrix and p i,u;j > 0 we conclude that for each (i, u) ∈ E × T ñ, λ i,u (t, η) > 0, iff t iu = 0.

(A16)

By relations (A14) and (A15) we have that

Lemma 4. 2 :

 2 Under condition (A2), i) the MLE η m (m; L) of the parameter η m , given by (34), is strongly consistent, that is, as L → ∞, we haveP α,η -a.s., that η m (m; L) → η m , ii) the sequence of random vectors ζ m (L) := √ L( η + m (m; L)η + m ), converges P α,η -weekly, as L → ∞, to a random vector ζ m , where ζ m ∼ N (0 , Γ

Corollary 4. 5 :

 5 Under condition (A2) and the probability measure P α,q , i) the MLE H(m; L) := H i (m; L) i∈E of the vector H(m) := H i (m) : i ∈ E is strongly consistent, and ii) the sequence of random vectors ζ m (L) := √ L( H(m; L) -H(m)) converges weakly, as L → ∞, to a random vector ζ m , where ζ m ∼ N (0 , Z (L) m ) and Z (L) m is a diagonal covariance matrix of dimension s × s given by

Yt

  n := 1 {Un-1≤ñ-1,Un=0} [t Zn-1Un-1Zn -E η (1 {Un=0} t Zn-1Un-1Zn | Z n-1 U n-1 )]. (A8)It is straightforward that the sequence of r.v. (Y n ) n∈N * is adapted to the filtration (F n ) n∈N * of the underlying probability space (Ω, A , P η ). Notice also that by constructionE η |Y n | < ∞ and E η (Y n | F n-1 ) = 0. Therefore, (Y n ) n∈N * is an (F n )-martingale difference. Since E η 1 {Un=0} t Zn-1Un-1Zn | Z n-1 U n-1 = j∈E-i p Zn-1,Un-1;j t Zn-1Un-1 j , iuj 1 {Zn-1=i,Un-1=u} (1 {Zn=j}p i,u;j ) = i,j∈E 2 * u∈Tñ t iuj N i,u (j; m) -N i,u (m)p i,u;j = tζ * (m). (A10)In order to apply the CLT for Martingales note thatE η (Y 2 n | F n-1 ) = 1 {Un-1≤ñ-1} Var η 1 {Un=0} t Zn-1Un-1Zn | Z n-1 U n-1 , 1≤ñ-1} Var η 1 {Un=0} t Zn-1Un-1Zn | Z n-1 , U n-1 = i∈E u∈Tñ N i,u (m)Var η 1 {Un=0} t iuZn | Z n-1 = i, U n-1 = u . (A11) Since E η 1 {Un=0} t 2 iuZn | Z n-1 = i, U n-1 = u = j∈E-i

  λ i,u (t, η)N i,u (m).

  {Un=0} t iuZn | Z n-1 = i, U n-1 = u =

	and
	E 2 η 1 j =i	t iuj p i,u;j	l =i	t iul p i,u;l ,
				=	t iuj p i,u;j	t iul p i,u;l ,	(A13)
					j =i		l =i
	by taking the difference of (A12) and (A13) to form the variance that appears in Relation
	(A11) we get
	m
	n=1	E η (Y 2 n | F n-1 ) =	i∈E u∈Tñ	N i,u (m)	j =i	l =i t iuj p i,u;j t iuj -	t iul p i,u;l
		=		N i,u (m)
			i∈E u∈Tñ		j,l∈E
					j,l =i
								t 2 iuj p i,u;j ,	(A12)

t iuj p i,u;j (δ jl t iujt iul p i,u;l ). (A14)

We denote by

λ i,u (t, η) := j,l∈E j,l =i t iuj p i,u;j (δ jl t iujt iul p i,u;l ). (

A15

)

Journal of Nonparametric Statistics

Acknowledgements

The authors are grateful to the editor, and two anonymous referees (including the associate editor). Their valuable comments and suggestions improved considerably this paper.

Appendix A. Proof of Results

Proof of Lemma 3.2: We define the random vector

where

Notice that if we denote by R(m) := diag{r iuj (m); i, j ∈ E, u ∈ T ñ}, where r iuj (m) = m N B i,u (m)

then it is easy to see that

Let also

Recall that the process N B i,u (m) is a counting process of the number of visits to the state (i, u) until time m of the irreducible and positive recurrent Markov chain (Z, U) (as a consequence of condition (A1)). From the ergodic theorem for positive recurrent MCs we have

Therefore by Relation (A5) and the expressions (A2) and (A4) we have

R(m)

Pη-a.s.

----→ R.

(A6)

If we establish that m -1/2 ζ(m) converges in distribution to a random vector ζ, then in virtue of the representation of ξ ñ(m) as in (A3) and the limiting result (A6), ξ ñ(m) will converge in distribution to the random vector Rζ. Note that by Relation (A1)

Since for (i, u) ∈ E × T ñ, we have π i,u > 0, then by (A16) we conclude that

Additionally, for arbitrary ǫ > 0

since Y m is bounded for any m ∈ N * . Since the limit expressed by equation ( A17) is strictly positive by (A18) and also (A19) holds true, we can apply the CLT for martingales to conclude by (A10) that P η -weakly,

By (A20), (A21) and the definition of ζ * we get that P η -weakly,

Relation (A22) is valid for any t ∈ R (s 2 -s)ñ , and consequently by (A3), (A6) and Slutsky theorem we get that P η -weakly,

The proof is complete for the minimal representation ξ * ñ(m), since the above covariance matrix coincides with Γ * ñ that results from Γ ñ which is given by Relation ( 22). In order to obtain ξ ñ(m) from the corresponding vector ξ * ñ(m) it is enough to add the components that correspond to the parameters p i,u;i , for all the couples (i, u). Since every such component is a linear combination of the others with the different j, identical to the case of the components of a multinomial distribution we conclude that the result can be extended easily in order to obtain the form of Γ ñ given by ( 22).

2 Proof of Proposition 3.3: We denote by Φ m the restriction of Φ, given by ( 12), to the domain that corresponds to the vector parameter η m in a minimal representation. Notice that

By Relation (A23) and delta method, since Φ ñ is an analytic function and therefore the total derivative Φ ′ ñ exists, we conclude that P η -weakly

The next step is to identify that the covariance matrix Φ ′ ñΓ * ñ(Φ ′ ñ) ⊤ coincides with ∆ ñ, which is given by ( 24). Firstly, it is easy to see that

where for i ∈ E,

Notice here that we consider the functions Φ ijk ordered with respect to the lexicographical order of (i, k, j). By Relation (12) we get

By combining (A24), ( A25) and (A26) we conclude that

We have that

Notice that

Hi( 1)

Hi( 1)

Hi( 2)

We have

and therefore

By equation ( A27) and (A28) we get for k < l + 1,

Since the matrices Φ ′ i;ñ Λ i (Φ ′ i;ñ ) ⊤ and K i are symmetric, their blocks coincide also for the case that k > l + 1. With similar arguments as in (A29), the diagonal blocks (k = l + 1) for the above matrices coincide.

2 Proof of the part iii) of Theorem 3.4: By (42), and since lim m→∞ m -1 π i,u|m = π i,u (ergodic theorem applied to the positive recurrent MC (Z, U)), we deduce that

and its inverse I -1 ∞ (η * ñ) equals Γ ñ in a minimal representation. ,u2;j2 (i1,u1,j1),(i2,u2,j2) , where the expectation is taken wrt to the stationary measure π of the MC (Z, U).

2 Proof of Corollary 4.5: Let R m : q + m → R m (q + m ) := H(m), where each component R i|m := H i (m) is given by (48). We conclude directly from the additive form of each component of R m that it is continuous and the total derivative exists. Therefore, since by the part i) of Theorem 4.4, the MLE q m (m; L) is strongly consistent, we can apply the continuous mapping theorem, and obtain the strong consistency of H(m; L). Likewise, by the part ii) of Theorem 4.4 and an application of delta method, the MLE H(m; L) is asymptotically normal with zero mean and diagonal covariance matrix of dimension s × s, given by

By ( 48) notice that R ′ i|m = 1 d * i|m , and consequently, by the form of K i given by ( 25), m given by (49).