Arithmetic Reduction of Age of order 1 ARD 1 Arithmetic Reduction of Degradation of order 1 ARD m

Arithmetic Reduction of Degradation of order m ARD ∞ Arithmetic Reduction of Degradation of order ∞ ARI 1 Arithmetic Reduction of Intensity of order 1 a.s. almost surely

Acronyms

Part I Introduction Safety and dependability are crucial issues in many industries (such as, e.g., railways, aircraft engines or nuclear power plants), which have lead to the development of the so-called reliability theory. For many years, only lifetime data were available and the rst reliability studies were focused on lifetime data analysis (see, e.g., [START_REF] Meeker | Statistical methods for reliability data[END_REF]), which still remains of interest in many cases. In that context and in case of repairable systems with instantaneous repairs, successive failure (or repair) times appear as the arrival points of a counting process, and failures hence correspond to recurrent events. As for the type of possible repairs, typical classical models are perfect (As-Good-As-New) and minimal (As-Bad-As-Old) repairs, leading to renewal and non homogeneous Poisson processes as underlying counting processes, respectively (see [START_REF] Barlow | Mathematical theory of reliability[END_REF]). Regarding the failure intensity, the eects of both these maintenance types are illustrated in Figure D.1: a minimal repair does not aect the failure intensity, while a perfect repair reduces it to its initial value. However, the reality often lies in-between, leading to the class of imperfect repairs (also represented in modeling, such as, e.g., virtual age models introduced by Kijima [START_REF] Kijima | Some results for repairable systems with general repair[END_REF], further studied in [START_REF] De Toledo | ARA and ARI imperfect repair models: Estimation, goodness-of-t and reliability prediction[END_REF]16], and extended in [START_REF] Brenière | Virtual age models with time-dependent covariates: A framework for simulation, parametric inference and quality of estimation[END_REF] where the authors introduce covariates in the virtual age model. Other possible models are geometric processes [START_REF] Lam | The geometric process and its applications[END_REF] (extended in [START_REF] Bordes | Extended geometric processes: semiparametric estimation and application to reliability[END_REF] and more recently in [START_REF] Dauxois | Semiparametric inference for an extended geometric failure rate reduction model[END_REF]) or, as already mentioned, models based on reduction of failure intensity [START_REF] De Toledo | ARA and ARI imperfect repair models: Estimation, goodness-of-t and reliability prediction[END_REF]16]. See, e.g., [START_REF] Doyen | On geometric reduction of age or intensity models for imperfect maintenance[END_REF] for a recent account and extensions of such models. See also [START_REF] Pham | Imperfect maintenance[END_REF] for more references and other models.

Nowadays, the development of online monitoring and the increasing use of sensors for safety assessment make it possible to get specic information on the health of a system and on its eective evolution over time, without waiting for the system failure. This information is often synthesized into a scalar indicator, which can for instance stand for the length of a crack, the thickness of a cable, the intensity of vibrations, corrosion level, ... This scalar indicator can be considered as a measurement of the deterioration level of the system. The evolution of this deterioration indicator over time is nowadays commonly modeled through a continuous-time and continuous-state stochastic process, which is often considered to have an increasing trend. Classical models include inverse Gaussian [START_REF] Ye | The Inverse Gaussian Process as a Degradation Model[END_REF] or Wiener processes (with trend) [START_REF] Hu | Optimum step-stress accelerated degradation test for wiener degradation process under constraints[END_REF][START_REF] Liu | A condition-based maintenance policy for degrading systems with age-and state-dependent operating cost[END_REF][START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF]. Also, the transformed Wiener process was lately introduced in [START_REF] Giorgio | A new age-and state-dependent degradation process with possibly negative increments[END_REF], and further studied in [START_REF] Giorgio | Bayesian estimation and prediction for the transformed Wiener degradation process[END_REF], where the degradation increments can be negative. All these stochastic processes are quite common in many other elds out of reliability theory, such as nance, insurance or epidemiology. This thesis focuses on gamma processes, which are widely used since they were introduced in the reliability eld mostly simultaneously by Abdel-Hameed [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF] and Çinlar [START_REF] Çinlar | Stochastic process for extrapolating concrete creep[END_REF]. This process is monotonous and is hence well adapted for modeling non decreasing degradation.

Before coming to the denition of a gamma process, let us start with the denition of the gamma distribution, which allows to set up the notations and parametrization used in this document.

A random variable X is said to be gamma distributed with shape parameter a > 0 and rate parameter b > 0 (denoted X ∼ Γ(a, b)), if its distribution admits the following probability density function (p.d.f.):

f X (x) = b a Γ(a)
x a-1 e -bx 1 R + (x) with respect to the Lebesgue measure. The corresponding expectation and variance are E(X) = a/b and V(X) = a/b 2 .

Let us now recall some well-known properties of the gamma distribution:

Let X 1 and X 2 be two independent and gamma distributed random variables with respective distributions Γ(a 1 , b) and Γ(a 2 , b), where a 1 , a 2 , b > 0. Then, for all c > 0, the random variables cX 1 and X 1 + X 2 are gamma distributed Γ (a 1 , b/c) and Γ(a 1 + a 2 , b) respectively. Now, let a(•) : R + → R + be a continuous and non decreasing function such that a(0) = 0 and let b > 0. Also let (X t ) t≥0 be a right-continuous stochastic process with left-side limits. Then, the stochastic process (X t ) t≥0 is a non homogeneous gamma process with shape function a(•) and rate parameter b, as soon as X 0 = 0 almost surely (a.s.); (X t ) t≥0 has independent increments; each increment is gamma distributed: for all 0 ≤ s < t, we have X t -X s ∼ Γ(a(t) -a(s), b), (see, e.g., [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF]). See [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF] and its references for a large overview of the gamma processes. Also, an example of simulated trajectories of a gamma process is given in Figure D. 2. In order to mitigate the degradation of the system over time and extends its lifetime, preventive maintenance actions can be considered, in addition to corrective repairs which are performed at failure.

In the context of deteriorating systems, many preventive maintenance policies from the literature consider condition-based maintenance (CBM) actions, where the preventive repair is triggered by the reaching of a preventive maintenance threshold by the deterioration level. In that context, "most of the existing CBM models have been limited to perfect maintenance actions", as noted by [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] (see also [START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF]). Some imperfect repair models are however emerging in the latest reliability literature, in this new context of deteriorating systems, see [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] for a recent review. Some models are based on the notion of virtual age previously introduced in the context of recurrent events (see, e.g., [START_REF] Giorgio | A new state-dependent degradation process and related model misidentication problems[END_REF][START_REF] Mercier | On the modelling of imperfect repairs for a continuously monitored gamma wear process through age reduction[END_REF]), where the system is rejuvenated by a maintenance action. Other models consider that an imperfect repair reduces the deterioration level of the system, such as [START_REF] Khatab | Condition-based selective maintenance for stochastically degrading multi-component systems under periodic inspection and imperfect maintenance[END_REF][START_REF] Letot | An adaptive degradation-based maintenance model taking into account both imperfect adjustments and agan replacements[END_REF][START_REF] Nicolai | Modelling and optimizing imperfect maintenance of coatings on steel structures[END_REF][START_REF] Ponchet | Maintenance policy on a nite time span for a gradually deteriorating system with imperfect improvements[END_REF][START_REF] Wu | A cost eective degradation-based maintenance strategy under imperfect repair[END_REF], which can be accompanied by some increase in the deterioration rate, as in [START_REF] Do | A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions[END_REF]. Also, some papers consider that the eciency decreases with the number of repairs (see, e.g., [START_REF] Liu | An imperfect maintenance policy for mission-oriented systems subject to degradation and external shocks[END_REF][START_REF] Zhang | An ameliorated improvement factor model for imperfect maintenance and its goodness of t[END_REF]), and further studies, as in [START_REF] Huynh | Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems[END_REF], deal with imperfect maintenance models such that (i) repairs have a random eciency (ii) the deterioration rate increases with the number of repairs.

In all these papers however, the main point mostly is on the optimization of a maintenance policy, including these imperfect maintenance actions together with perfect repairs (replacements). Up to our knowledge, very few papers from the literature deal with statistical issues concerning imperfect repair models for deteriorating systems, except from [START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF], where the authors suggest a maximum likelihood method for estimating the parameters of the Wiener process (that describes the deterioration out of repairs) together with an iterative procedure based on a Kalman lter for the dierent factors implied in successive imperfect repairs.

In the context of systems subject to deterioration and imperfect repairs, the estimates for the parameters of the underlying degradation process and maintenance eciency are of great use for maintenance policies optimization. Indeed, once the parameters have been estimated, the future behavior of the maintained system can be predicted, which allows to adapt (optimize) the periodicity of the maintenance actions and eciently plan a general overhaul for instance. From a safety point of view, the principal inquiry is to ensure that the maintenance actions are eective enough to keep with a high probability the degradation level below a xed threshold (safety level). As long as this safety level is not reached, the maintenance actions can be adjusted, either by adapting their periodicity or by improving their eciency (if possible). Of course, apart from the previous safety concern, the maintenance costs are another issue.

As an example, in [START_REF] Wu | A cost eective degradation-based maintenance strategy under imperfect repair[END_REF], the costs minimization is based on the monitoring time and on the imperfect maintenance eciency. In [START_REF] Huynh | Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems[END_REF], the author considers a threshold for the degradation level, beyond which an imperfect maintenance action is performed. The optimization is made with respect to this threshold and the inspections periodicity. Finally, in [START_REF] Wu | Maintenance policy for a system with a weighted linear combination of degradation processes[END_REF], a maintenance policy is proposed, where a replacement is performed either when the degradation exceeds a given threshold or when a xed number of imperfect preventive maintenance actions are conducted. See, e.g., the three papers cited above and their reference for an overview on maintenance policies optimization.

This thesis focuses on the development and applications of estimation procedures for three specic imperfect repair models in the context of a gamma deteriorating system. The document is split into four parts, including the present introduction.

Part II deals with the Arithmetic Reduction of Degradation models of order one and innity (ARD 1 and ARD ∞ ), where each maintenance action reduces the deterioration level of the system. The ARD 1 model was rst introduced in [START_REF] Castro | Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs[END_REF] and further studied in [34]. Mimicking the Arithmetic Reduction of Intensity (ARI 1 ) model of order 1 developed by [16] in the context of recurrent events, the idea of this model is that a maintenance action removes a proportion ρ of the degradation accumulated by the system from the last maintenance action (where ρ ∈ [0, 1)). Based on the same idea, [16] also dened the ARI model of innite order for recurrent events, that we extend here to the degradation framework in order to introduce the ARD ∞ model. Regarding this imperfect repair model, each maintenance action removes a proportion ρ of the current degradation, that is the degradation accumulated by the system from the initial time t = 0. Once these models dened, we place ourselves in a fully parametric framework and the observation scheme is stated. The Moments Estimation (ME) and Maximum Likelihood Estimation (MLE) methods are developed in Chapter 2 for both models. To be more precise, the identiability is studied and then we provide expressions for the parameters estimates. Numerical experiments based on simulated data are conducted in Chapter 3. In Chapter 4, we propose an original estimator for ρ, which does not depend on the underlying gamma process parameters, leading to a semiparametric framework. The idea of this estimator has come from a preliminary study in the framework based on the MLE method, where we observed that for one single trajectory, the minimum of the set of admissible ρ's has quite an interesting behavior when the shape function a of the underlying Gamma deterioration process is concave, getting quickly very close to the unknown eciency parameter when the number of repairs increases. This semiparametric estimate was rst developed in the context of the ARD 1 model, leading to a publication (see [38]) which is exactly reproduced here in Sections 4.1 to 4.6. In Section 4.7 we propose an extension of this work to the ARD ∞ model, which is specic to this thesis.

In Part III, we consider an imperfect repair model based on the virtual age introduced by [34]: the Arithmetic Reduction of Age model of order one (ARA 1 ). Following the idea of [16] in the context of lifetime analysis, this model is such that a maintenance action removes a proportion ρ of the age of the system accumulated since the last maintenance action. Unlike for the ARD 1 model, here the system is rejuvenated, that is it is put back to the exact situation where it was some time before. The rst steps of the work are similar to those of the previous part: the model as well as the observation scheme are dened, the ME method is studied from the point of view of the model parameters identiability, and estimates expressions are provided. Then, the MLE method is developed. However, due to dependency issues, the likelihood function appears to be a product of integrals of large dimension, and thus numerical estimations becomes dicult to compute in a classical way. Hence, the MLE method requires approximating the integrals by the Monte Carlo and randomized Quasi Monte Carlo methods. Also, in order to avoid the computation of high-dimensional integrals, alternatives to the maximum likelihood method are developed: the composite maximum likelihood and the half data methods. The rst one is based on the composite marginal likelihood constructed under independence assumptions, also referred to in the literature as the independence likelihood (see, e.g., [START_REF] Chandler | Inference for clustered data using the independence loglikelihood[END_REF]). In the present work, this method consists in considering that the observations are independent and thus obtain an approximation of the likelihood function, which allows to develop the maximum composite likelihood estimation method. See [START_REF] Varin | Composite likelihood methods[END_REF] for an overview of composite likelihood methods. The half data method is based on sub-sample methods introduced by [START_REF] Hartigan | Using subsample values as typical values[END_REF] and commonly used within bootstrap estimation (see, e.g., [START_REF] Babu | Subsample and half-sample methods[END_REF]). Here only one out of two observations are taken into account, which eliminates the dependence issues. Once again, the study of the estimation methods is based on the parameters identiability, which was very challenging in one case, as well as the expressions of either the estimates or the log-likelihood. However, the identiability is not treated regarding two out of the six methods we studied since the simulation studies reveal identiability issues for these specic methods. Lastly in Chapter 7, some illustrations of the numerical performances of these methods are provided in two stages. A rst study allows to eliminate methods with poor performance while a second study allows to select, among the remaining methods, the most appropriate one as a function of the observations characteristics.

Finally, Part IV highlights the conclusions of the thesis and future perspectives.

Part II

Imperfect repairs based on reduction of the degradation level

Chapter 1

Introduction 1.1 Preliminary

A system is considered whose intrinsic deterioration is modelled by a gamma process (X t ) t≥0 with shape function a(•) and scale parameter b, as dened in the introduction. The system is subject to periodic (period T ) and instantaneous imperfect maintenance actions, where each maintenance action removes ρ% of the deterioration accumulated from the last m maintenance actions, where m ∈ N * and ρ in [0, [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF] are xed. This model is called Arithmetic Reduction of Degradation model of order m (ARD m ). The Euclidean parameter ρ is a measure of the eciency of the imperfect maintenance action.

In the following, two models are considered: the ARD 1 model and the ARD ∞ model. The denitions of both models are rst given in the two following sections. Two parametric estimation methods are developed in Chapter 2, and numerical experiments based on simulated data are provided in Chapter 3. Finally, a semiparametric estimate of the maintenance actions eciency is developed in Chapter 4.

Arithmetic Reduction of Degradation model of order one

Regarding this model, a maintenance action reduces the degradation accumulated from the last maintenance action only. The repairs have an eect (or eciency) dened by a parameter ρ ∈ [ 0 , 1) and thus, at time jT , j ∈ N * , the maintenance eect results in a reduction of a proportion ρ of the deterioration accumulated over the time interval [ (j -1)T , jT ).

Let (Y t ) t≥0 stands for the random process describing the degradation level evolution of the maintained system according to the ARD 1 model. Let (X (j) ) j∈N * be a sequence of independent copies of (X t ) t≥0 , where X (j) describes the intrinsic deterioration of the system over ((j -1)T, jT ). The system is assumed to be in the perfect working order at the initial time, that is Y 0 = X

(1) 0 = 0. Over the time interval [0 , T ), the system deteriorates according to X

t . Therefore:

Y t = X (1) t and Y T = (1 -ρ) X (1)
T -X

(1) 0

= (1 -ρ)X (1)
T .

For t ∈ [T , 2T ), the system's degradation level is equal to the sum of the degradation accumulated between times t and T , that is

X (2) t -X (2)
T , and the degradation level right after the rst maintenance Chapter 1. Introduction action, thus:

Y t = Y T + X (2) t -X (2)
T

and

Y 2T = Y 2T -+ ρ X (2)
2T --X

(2)

T = (1 -ρ) X (2)
2T -X

(2) T + X

T -X

(1) 0

.

Hence, for any t in [ nT , (n + 1)T ), with n in N * , we have:

Y nT = (1 -ρ) n j=1 X (j)
jT -X (j) (j-1)T

and

Y t = Y nT + X (n+1) t -X (n+1) nT .
Note that the random variables Y nT and X (n+1) t -X (n+1) nT are independent and gamma distributed, Γ (a(nT ) , b/(1 -ρ)) and Γ (a(t) -a(nT ) , b), respectively. This denition as well as the following property are already given in [34], and we recall it here for sake of completeness.

Proposition 1. For each n in N and nT ≤ t < (n + 1)T , the expectation of the degradation level Y t is given by E (Y t ) = (a(t) -ρa(nT ))/b and its variance by V (Y t ) = (a(t) -ρ (2 -ρ) a(nT ))/b 2 .

Proof. Let us set n in N and nT ≤ t < (n + 1)T . The expression of Y t is the following

Y t = Y nT + X (n+1) t -X (n+1) nT ,
where Y nT and X (n+1) t -X (n+1) nT are independent and gamma distributed Γ (a(nT ) , b/(1 -ρ)) and Γ (a(t) -a(nT ) , b) respectively. Hence, the expectation as well as the variance of Y t can be expressed as the sum of both terms, and thus

   E(Y t ) = E(Y nT ) + E X (n+1) t -X (n+1) nT V(Y t ) = V(Y nT ) + V X (n+1) t -X (n+1) nT
We now replace the expressions of the mean and variance of the gamma distribution, see Part I, which leads to

     E(Y t ) = (1 -ρ)a(nT ) + a(t) -a(nT ) b V(Y t ) = (1 -ρ) 2 a(nT ) b 2 + a(t) -a(nT ) b 2
and nishes the proof. represents the degradation evolution expected of the maintained system, as well as that of the intrinsic deterioration expected (unmaintained system).

Arithmetic Reduction of Degradation model of order innity

This model is very close to the previous one. The only dierence is that a maintenance action aect the degradation accumulated from the initial time t = 0. Thus at times jT , for j in N * , the maintenance action results in a reduction of a proportion ρ of the current degradation level.

Let (Y t ) t≥0 be the random process describing the degradation level evolution of a system according to the ARD ∞ model. Let (X (j) ) j∈N * be a sequence of independent copies of (X t ) t≥0 . We consider once again that Y 0 = X

(1) 0 = 0. Both ARD 1 and ARD ∞ models have identical behaviour over the time interval [0 , 2T ). Then, for all t in [0 , T ) we have:

Y t = X (1) t and Y T = (1 -ρ) X (1) T -X (1) 0 = (1 -ρ)X (1) T ,
and for all t in [ T , 2T )

Y t = Y T + X (2) t -X (2) T .
When the second maintenance action happens, a proportion ρ of the degradation accumulated from the

time t = 0 is removed, thus Y 2T = (1 -ρ) Y T + X (2) 2T -X (2) T
, which can be written as:

Y 2T = (1 -ρ) 2 X (1) T -X (1) 0 + (1 -ρ) X (2) 2T -X (2) T .
From this, for any t in [nT , (n + 1)T ) with n in N, we easily derive the expression of Y t , which is: in N * and 1 ≤ i ≤ n, we have:

Y t = n j=1 (1 -ρ) n-j+1 X (j) jT -X (j) (j-1)T + X (n+1) t -X (n+1) nT
X (n+1) t -X (n+1) nT ∼ Γ (a(t) -a(nT ) , b) ; (1 -ρ) n-j+1 X (j) jT -X (j) (j-1)T ∼ Γ a(jT ) -a((j -1)T ) , b (1-ρ) n-j+1 .
Hence, the mean and variance of Y t are given by

E (Y t ) = 1 b   n j=1 (1 -ρ) n-j+1 (a(jT ) -a((j -1)T )) + (a(t) -a(nT ))   and V(Y t ) = 1 b 2   n j=1 (1 -ρ) 2(n-j+1) (a(jT ) -a((j -1)T )) + (a(t) -a(nT ))   for all t ∈ [nT , (n + 1)T ) with n ∈ N.
Based on the same framework as in Example 1, Figure 1.3 represents the degradation evolution expected of the maintained system, as well as that of the intrinsic deterioration mean (unmaintained system). Figure 1.4 represents the variances of the same quantities.

For sake of simplicity, from now on, we set ∆X

(j) = X (j) jT -X (j)
(j-1)T , which represents the increment of the intrinsic deterioration over [(j -1)T, jT ).

Chapter 2

Parametric inference for the Arithmetic Reduction of Degradation models

Preliminary

In this chapter, the Moments Estimation (ME) and the Maximum Likelihood Estimation (MLE) methods are developed in the framework of both ARD models dened above, in order to estimate the models parameters. The periodicity T is assumed to be known, thus the parameters of interest are the parameters of the shape function, the scale parameter b and the maintenance eciency ρ. The estimation methods are then tested in the next chapter, within the framework of a power law shape function. To be more precise, the shape function is dened as a : t → αt β with α, β > 0. We set θ = (α, β, ρ, b) ∈ Θ the parameter set, with Θ = (0, ∞) 2 × [0, 1), ×(0, ∞).

We assume that the degradation level is measured right before the rst n maintenance actions for some n in N * , that is at times T -, 2T -,. . . ,nT -.. Also, s i.i.d. systems are considered. For sake of readability, let us dene the following notations:

Y j = Y jT -for 1 ≤ j ≤ n; Y = Y jT -1≤j≤n = (Y j ) 1≤j≤n ; y (i) j
is the observed degradation level of the ith maintained system at times jT -for 1 ≤ j ≤ n and 1 ≤ i ≤ s, which is a realisation of the r.v. Y j ;

y (i) = y (i) j 1≤j≤n
is the complete observations set related to the ith maintained system, it is a realisation of Y ;

y is the complete observations set, that is y = y (i) j 1≤j≤n, 1≤i≤s
, it corresponds to s i.i.d. realisations of Y.

Moments method estimation 2.2.1 Description of the method

In order to estimate the parameter θ by the ME method, the distance between the empirical moments and the theoretical moments must be minimized. Here the centered moments are considered, except for Chapter 2. Parametric inference for the Arithmetic Reduction of Degradation models the rst moment. Let us consider the distance function D which is dened as

D(θ , θ 0 ) = d k=1 n j=1 m k (θ , jT -) -m k (θ 0 , jT -) 2
where θ 0 is the true value of the parameter, d the maximal order of the moments, n the number of observations by trajectory and m k (θ , jT -) the kth moment (centered if k > 1) at time jT -, that is

m 1 (θ, jT -) = E Y jT -and ∀k ≥ 2, m k (θ, jT -) = E Y jT --E Y jT - k .
Then, an estimation θ of θ 0 can be obtained through

θ = arg min θ∈Θ D(θ)
with D(θ) the empirical version of D(θ , θ 0 ), where true unknown moments are replaced by their empirical counterparts. This empirical version of D has the following expression

D(θ) = d k=1 n j=1 m k (θ , jT -) -mk (jT -) 2 (2.1)
with mk (jT -) the kth empirical moment (centered if k > 1) at time jT -. These empirical moments are empirical estimates of the moments m k (θ 0 , jT -) and are dened by

m1 (jT -) = 1 s s i=1 y (i) j and ∀k ≥ 2, mk (jT -) = 1 s s i=1 y (i) j -m1 (jT -) k .
This method is based on the classical approach of the Generalized Method of Moments as exposed in [START_REF] Liu | Generalized method of moment estimation of multivariate multifractal models[END_REF]. The two following sections are devoted to the application of the ME method as dened above, in the case of the ARD 1 and ARD ∞ models, beginning by the model of order 1.

Application to the Arithmetic Reduction of Degradation model of order one

Let us recall that the shape function is dened by a(t) = αt β with α, β > 0, and thus the estimation

focuses on four parameters θ = (α , β , ρ, b) ∈ Θ with Θ = (0 , +∞) 2 × [0 , 1) × (0 , +∞).
In order to apply the ME method, we rst need to see how many moments d and how many data n are necessary to identify the model parameters, which is done in next proposition.

Proposition 2. The parameters of the ARD 1 model are identiable from the ME method, that is D(θ , θ 0 ) = 0 implies that θ = θ 0 for all θ, θ 0 ∈ Θ, as soon as d ≥ 2 and n ≥ 3. In other words, the identiability holds from the ME method as soon as the systems are observed at times T -, 2T -and 3T -, and if at least the rst two moments (expectation and variance) are used.

Proof. Here we prove that the assertion D(θ , θ 0 ) = 0 ⇒ θ = θ 0 is true ∀ θ, θ 0 ∈ Θ as soon as n ≥ 3 and d ≥ 2. Assume that D(θ , θ 0 ) = 0, and let θ = (α , β , b , ρ, b) and θ 0 = (α 0 , β 0 , ρ 0 , b 0 ) be in Θ. Based on the denition of the function D, we have for all j in {1 , ... , n}:

m 1 (θ , jT -) = m 1 (θ 0 , jT -) m 2 (θ , jT -) = m 2 (θ 0 , jT -)
When j = 1, the rst two moments are those of the random variable X

T , hence this system is equivalent to

α b T β = α0 b0 T β0 α b 2 T β = α0 b 2 0 T β0 thus b = b 0 αT β = α 0 T β0
Note that only the parameter b is identied. Taking into account the second observation time, that is when j = 1, 2 we can write

           α b T β = α0 b0 T β0 α b 2 T β = α0 b 2 0 T β0 α b (2T ) β -α b ρT β = α0 b0 (2T ) β0 -α0 b0 ρ 0 T β0 α b 2 (2T ) β -α b 2 ρ(2 -ρ)T β = α0 b 2 0 (2T ) β0 -α0 b 2 0 ρ 0 (2 -ρ 0 )T β0
and thus

           α b T β = α0 b0 T β0 α b 2 T β = α0 b 2 0 T β0 α b T β 2 β -ρ = α0 b0 T β0 2 β0 -ρ 0 α b 2 T β 2 β -ρ(2 -ρ) = α0 b 2 0 T β0 2 β0 -ρ 0 (2 -ρ 0 )
which can be reduced to

           b = b 0 αT β = α 0 T β0 2 β -ρ = 2 β0 -ρ 0 2 β -ρ(2 -ρ) = 2 β0 -ρ 0 (2 -ρ 0 )
Now let us rst assume that ρ 0 = 1 -ρ. The third and fourth equations lead to ρ(1

-ρ) = ρ 0 (1 -ρ 0 ),
hence we have ρ = ρ 0 . This implies that 2 β = 2 β0 and thus β = β 0 from the third equation. Finally, the second equation provides α = α 0 . Now let us assume that ρ 0 = 1 -ρ. In that case, we need to add the third observation, leading to the following expressions for the expectation and variance when j = 1, 3:

           α b T β = α0 b0 T β0 α b 2 T β = α0 b 2 0 T β0 α b T β 3 β -2 β ρ = α0 b0 T β0 3 β0 -2 β0 ρ 0 α b 2 T β 3 β -2 β ρ(2 -ρ) = α0 b 2 0 T β0 3 β0 -2 β0 ρ 0 (2 -ρ 0 )
This system can be written as Chapter 2. Parametric inference for the Arithmetic Reduction of Degradation models

           b = b 0 αT β = α 0 T β0 3 β -2 β ρ = 3 β0 -2 β0 ρ 0 3 β -2 β ρ(2 -ρ) = 3 β0 -2 β0 ρ 0 (2 -ρ 0 )
and because ρ 0 = 1 -ρ by assumption, we have

           b = b 0 αT β = α 0 T β0 3 β -2 β ρ = 3 β0 -2 β0 + 2 β0 ρ 3 β -2 β ρ(2 -ρ) = 3 β0 -2 β0 + 2 β0 ρ 2
The third and fourth equations lead to

2 β (ρ -ρ 2 ) = 2 β0 (ρ -ρ 2 ). If ρ = 0 then ρ 0 = 1 -ρ = 1
which is impossible since ρ 0 ∈ [0, 1). Hence ρ ∈ (0, 1), which implies that 2 β = 2 β0 and thus β = β 0 , and nally α = α 0 from the second equation. As a result, from the third equation we deduce that ρ = ρ 0 .

Note that because ρ 0 = 1 -ρ, necessarily ρ = ρ 0 = 0.5.

Hence, only the rst three observations and the rst two moments are necessary in order to obtain the identiability of the model parameters.

In the following, we consider d = 2, thus Equation (2.1) becomes

D(θ) = n j=1 m 1 (θ , jT -) -m1 (jT -) 2 + m 2 (θ , jT -) -m2 (jT -) 2
which can be written as

D(θ) = n j=1 E Y jT --m1 (jT -) 2 + V Y jT --m2 (jT -) 2 .
From now on, another parametrization of the empirical distance function D(θ) is considered in order to simplify our calculations. Let us set θ = (µ , η , b , ρ), with µ = α/b and η = α/b 2 . Thus we have for all j in {1 , ... , n},

E(Y jT -) = µT β j β -ρ (j -1) β and V(Y jT -) = ηT β j β -ρ(2 -ρ) (j -1) β , which implies
D(θ ) = µ 2 g 1 (β , ρ) -2µ g 2 (β , ρ) + η 2 h 1 (β , ρ) -2η h 2 (β , ρ) + C (2.2)
where the functions g 1 , g 2 , h 1 and h 2 are independent of µ and η, and C is a constant independent of θ .

To be more precise, the quantity C is given by

C = n j=1 m1 (jT -) 2 + m2 (jT -) 2
2.2. Moments method estimation and the functions g 1 , g 2 , h 1 and h 2 have the following expressions

g 1 (β , ρ) = T 2β   n j=1 j 2β -2ρ n j=1 (j (j -1)) β + ρ 2 n j=1 (j -1) 2β   g 2 (β , ρ) = T β   n j=1 j β m1 (jT -) -ρ n j=1 (j -1) β m1 (jT -)   h 1 (β , ρ) = T 2β   n j=1 j 2β -2ρ(2 -ρ) n j=1 (j (j -1)) β + ρ 2 (2 -ρ) 2 n j=1 (j -1) 2β   h 2 (β , ρ) = T β   n j=1 j β m2 (jT -) -ρ(2 -ρ) n j=1 (j -1) β m2 (jT -)  
We now look for critical points of D(θ ), by looking for the zeros of the partial derivatives of D with respect to θ. However, due to the expression of D(θ ) given by Equation (2.2), we only write down

∂ µ D(θ ) and ∂ η D(θ )
as the other partial derivatives are of no use, because of their complexity. This leads to:

∂ µ D(θ ) = 0 ∂ η D(θ ) = 0 which provides 2µ g 1 (β , ρ) -2g 2 (β , ρ) = 0 2η h 1 (β , ρ) -2h 2 (β , ρ) = 0
Hence we can express the parameters µ and η with respect to β and ρ as follows:

µ = h µ (β, ρ) = g 2 (β , ρ) g 1 (β , ρ) and η = h η (β, ρ) = h 2 (β , ρ) h 1 (β , ρ)
. Based on this, Equation (2.2) can be now written as

D(β, ρ) ≡ D (h µ (β, ρ), h η (β, ρ), β, ρ) = C - g 2 2 (β , ρ) g 1 (β , ρ) + h 2 2 (β , ρ) h 1 (β , ρ) .
Then, the function D(β, ρ) is numerically minimized to obtain an estimation of β and ρ. Thereafter we estimate µ and η by plugging, setting μ = h µ ( β , ρ) and η = h η ( β , ρ). From this we derive expressions of α and b, and nally we obtain the following estimate

θ = g 2 ( β , ρ)h 1 ( β , ρ) g 1 ( β , ρ)h 2 ( β , ρ) , β , ρ, g 2 2 ( β , ρ)h 1 ( β , ρ) g 2 1 ( β , ρ)h 2 ( β , ρ) . (2.3)
where

β , ρ = arg min (β,ρ)∈(0,∞)×[0,1)

D(β, ρ)

We now look at the ME method in the case of an ARD ∞ model.

Application to the Arithmetic Reduction of Degradation model of order innity

This case is similar to the previous one, thus less details are given in this section, beginning with the following proposition about the parameters identiability.

Proposition 3. The parameters of the ARD 1 model are identiable from the ME method, that is D(θ , θ 0 ) = 0 implies that θ = θ 0 for all θ, θ 0 ∈ Θ, as soon as d ≥ 2 and n ≥ 3.

Proof. Here we prove that the assertion D(θ , θ 0 ) = 0 ⇒ θ = θ 0 is true ∀ θ, θ 0 ∈ Θ as soon as n ≥ 3 and d ≥ 2. Assume that D(θ , θ 0 ) = 0, and let θ = (α , β , b , ρ, b) and θ 0 = (α 0 , β 0 , ρ 0 , b 0 ) be in Θ. So given the function D, we have for all j in {1 , ... , n}:

m 1 (θ , jT -) = m 1 (θ 0 , jT -) m 2 (θ , jT -) = m 2 (θ 0 , jT -)
Note that the mean and expectation of the degradation level Y t over [0, 2T ) is equivalent for both ARD 1

and ARD ∞ model. Hence, most of the present proof is similar to this of Proposition 2, hence it can be concluded that the identiability holds if n = 2 and ρ = 1 -ρ 0 . Now assume that ρ 0 = 1 -ρ and let us add the third observation (j = 3). Let us recall that for j = 1,

we have b = b 0 and αT β = α 0 T β0 . Hence, we derive the following equations when j = 1, 2, 3:

                     b = b 0 αT β = α 0 T β0 2 β -ρ = 2 β0 -ρ 0 2 β -ρ(2 -ρ) = 2 β0 -ρ 0 (2 -ρ 0 ) 3 β -2 β + (1 -ρ)(2 β -1) + (1 -ρ) 2 = 3 β0 -2 β0 + (1 -ρ 0 )(2 β0 -1) + (1 -ρ 0 ) 2 3 β -2 β + (1 -ρ) 2 (2 β -1) + (1 -ρ) 4 = 3 β0 -2 β0 + (1 -ρ 0 ) 2 (2 β0 -1) + (1 -ρ 0 ) 4
We have ρ 0 = 1 -ρ by assumption, thus

                     b = b 0 αT β = α 0 T β0 2 β -ρ = 2 β0 -(1 -ρ) 2 β -ρ(2 -ρ) = 2 β0 -(1 -ρ)(1 + ρ) 3 β -2 β + (1 -ρ)(2 β -1) + (1 -ρ) 2 = 3 β0 -2 β0 + ρ(2 β0 -1) + ρ 2 3 β -2 β + (1 -ρ) 2 (2 β -1) + (1 -ρ) 4 = 3 β0 -2 β0 + ρ 2 (2 β0 -1) + ρ 4
The third and fourth equations are identical and provides

2ρ -1 = 2 β -2 β0 (2.4)
while the last two ones leads to

ρ(1 -ρ) ρ 2 -3ρ + 2 β + 1 = ρ(1 -ρ) ρ 2 + ρ + 2 β0 -1 .
Note that if ρ = 0, then ρ 0 = 1 which is impossible since ρ 0 belongs to [0, 1). Hence ρ ∈ (0, 1), and the 2.2. Moments method estimation equation above is equivalent to

4ρ -2 = 2 β -2 β0 . (2.5)
Finally, from Equations (2.4) and (2.5), we derive that 4ρ -2 = 2ρ -1. Therefore ρ = 1/2 and ρ 0 = 1 -ρ = 1/2, and nally ρ = ρ 0 . This leads to β = β 0 and then α = α 0 , which nishes the proof.

The function D (see Equation (2.1)) is also considered here, and an estimation of θ will be obtained through the minimization of this function. The same parametrisation θ = (µ , η , b , ρ) is used, with µ = α/b and η = α/b 2 , and for j in {1 , ... , n}, we have:

E(Y jT -) = µg (j) 1 (β , ρ) and V(Y jT -) = ηg (j) 2 (β , ρ), with g (j) 1 (β , ρ) = T β j β -(j -1) β + j-1 i=1 (1 -ρ) j-i i β -(i -1) β g (j) 2 (β , ρ) = T β j β -(j -1) β + j-1 i=1 (1 -ρ) 2(j-i) i β -(i -1) β .
Thus the function D(θ ) can be written as

D(θ ) = n j=1 µg (j) 1 (β , ρ) -m1 (jT -) 2 + n j=1 ηg (j) 2 (β , ρ) -m2 (jT -) 2 .
(2.6)

Once again we look for the zeros of the partial derivatives of D with respect to µ and η. Here we have

∂ µ D(θ ) = 0 ∂ η D(θ ) = 0 which is equivalent to            µ n j=1 g (j) 1 (β , ρ) 2 - n j=1 g (j) 1 (β , ρ) m1 (jT -) = 0 η n j=1 g (j) 2 (β , ρ) 2 - n j=1 g (j) 2 (β , ρ) m2 (jT -) = 0
Based on the above equations we can express µ and η as functions of β and ρ that cancel the partial derivatives of D with respect to µ and η, leading to

µ = h µ (β , ρ) = n j=1 g (j) 1 (β , ρ) m1 (jT -) n j=1 g (j) 1 (β , ρ) 2
Chapter 2. Parametric inference for the Arithmetic Reduction of Degradation models and

η = h η (β , ρ) = n j=1 g (j) 2 (β , ρ) m2 (jT -) n j=1 g (j) 2 (β , ρ) 2 .
Then, replacing µ and η by h µ (β , ρ) and h η (β , ρ) in Equation (2.6), we obtain a function D which only depends on β and ρ, and is given by

D(β, ρ) ≡ D (µ(β, ρ), η(β, ρ), β, ρ) = C -            n j=1 g (j) 1 (β , ρ) m1 (jT -)   2 n j=1 g (j) 1 (β , ρ) 2 +   n j=1 g (j) 2 (β , ρ) m2 (jT -)   2 n j=1 g (j) 2 (β , ρ) 2         
where C is the same constant as dened in the previous section. Finally, the minimization of this function with respect to β and ρ provides estimations of β and ρ, and thus of the entire parameter set θ . Then, from the relationship between θ and θ we derive an estimate of θ, which is given by

θ =          n j=1 g (j) 1 ( β , ρ) m1 (jT -) n j=1 g (j) 2 ( β , ρ) 2 n j=1 g (j) 1 ( β , ρ) 2 n j=1 g (j) 2 ( β , ρ) m2 (jT -) , β , ρ,   n j=1 g (j) 1 ( β , ρ) m1 (jT -)   2 n j=1 g (j) 2 ( β , ρ) 2   n j=1 g (j) 1 ( β , ρ) 2   2 n j=1 g (j) 2 ( β , ρ) m2 (jT -)          .

Maximum likelihood estimation 2.3.1 Application to the Arithmetic Reduction of Degradation model of order one

In this section we deal with the MLE method. We begin with the case where one single system is observed (s = 1) and we start by computing the joint p.d.f. of the random vector Y , from where the likelihood function is easily expressed for s > 1. With that aim, we need the following technical result.

Proposition 4. The degradation level Y jT right after the maintenance at time jT can be expressed with respect to the random variables Y 1 , . . . , Y j and the maintenance eciency parameter ρ as follows:

Y jT = (1 -ρ) j k=1 ρ j-k Y k , for all 1 ≤ j ≤ n.

Maximum likelihood estimation

Proof. The result is proved by induction on j. When j = 1, the result is veried as (1 -ρ)

j k=1 ρ j-k Y k = (1 -ρ)Y 1 = Y T . Now assume that Y jT = (1 -ρ) j k=1 ρ j-k Y k , for some 1 ≤ j ≤ n,
we want to prove that

Y (j+1)T = (1 -ρ) j+1 k=1 ρ j-k+1 Y k ,
or equivalently that j+1) , which can also be expressed as

Y (j+1)T = (1 -ρ)Y j+1 + (1 -ρ) j k=1 ρ j-k+1 Y k . (2.7) Recall that Y (j+1)T = Y jT + (1 -ρ)∆X (
Y (j+1)T = Y jT + (1 -ρ) (Y j+1 -Y jT ) .
This leads to

Y (j+1)T = (1 -ρ)Y j+1 + ρY jT .
By assumption on Y jT , we have

Y (j+1)T = (1 -ρ)Y j+1 + ρ(1 -ρ) j k=1 ρ j-k Y k
which can be written as Equation (2.7) and thus nishes the proof.

Lemma 1. Let y = (y 1 , . . . , y n ) be a realisation of Y . The p.d.f. f Y of Y is given by

f Y (y) = n j=1 f ∆X (j) (g j (ρ, y)) ,
where f ∆X (j) is the p.d.f. of the random variable ∆X (j) , and f Yj |(Y1, ... ,Yj-1) (y j | (y 1 , . . . , y j-1 )) .

       g j (ρ, y) = y j -(1 -ρ) j-1 k=1 ρ j-k-1 y k g 1 (ρ, y) = y 1
Then, we have by denition Y j f Yj |(Y1, ... ,Yj-1) (y j | (y 1 , . . . , y j-1 )) = f Y (j-1)T +∆X (j) |(Y1, ... ,Yj-1) (y j | (y 1 , . . . , y j-1 )) .

Let us now recall that for all 2 ≤ j ≤ n we have ∆X (j) ⊥ (Y 1 , . . . , Y j-1 ) and

Y (j-1)T = (1 -ρ) j-1 k=1 ρ j-k-1 Y k from Proposition 4. We derive that, f Y (j-1)T +∆X (j) |(Y1, ... ,Yj-1) (y j | (y 1 , . . . , y j-1 )) = f ∆X (j) (g j (ρ, y)) .
This nishes the proof as f ∆X (1) (g 1 (ρ, y)) = f Y1 (y 1 ) by denition of Y 1 and g 1 .

We now consider s independent and identical systems. The previous lemma allows us to write down the expression of the likelihood function as follows:

L (θ|y) = s i=1 f Y (y (i) ) (2.8)
which is equivalent to

L (θ|y) = s i=1 n j=1 f ∆X (j) g (i) j (ρ) (2.9) with        g (i) j (ρ) = y (i) j -(1 -ρ) j k=1 ρ j-k-1 y (i) k g (i) 1 (ρ) = y (i) 1
Let us recall that for 1 ≤ j ≤ n, the random variables ∆X (j) are gamma distributed with shape and scale parameters ∆A (j) and b, respectively. As a result, Equation (2.9) becomes

L (θ|y) = s i=1 n j=1
b ∆A (j) Γ ∆A (j) g

(i) j (ρ) ∆A (j) -1 e -bg (i) j (ρ) (2.10)
where ∆A (j) = αT β j β -(j -1) β .

Maximum likelihood estimation

For similar reasons as for parameters identiability when applying the ME method, we now study the identiability in the case of the MLE method.

Proposition 5. The parameters of the ARD 1 model are identiable from the likelihood function, that is L (θ|y) = L (θ 0 |y) for all y = (y 1 , . . . , y n ) implies that θ 0 = θ and θ, θ 0 ∈ Θ, as soon as n ≥ 2. In other words, the identiability holds from the likelihood function as soon as observations are conducted twice, at times T -and 2T -, and whatever s is.

Proof. Only the case where one single system is observed (s = 1), which is sucient to conclude. Assume that L (θ|y) = L (θ 0 |y) for all y = (y 1 , . . . , y n ), the purpose is to show that θ = θ 0 for all θ = (α, β, ρ, b)

and θ 0 = (α 0 , β 0 , ρ 0 , b 0 ) in Θ as soon as n ≥ 2, whatever s is.

First, when n = 1, the observation is a realization of ∆X (1) , which does not depend on rho, so that the model obviously is non identiable. Assume that n = 2. We have by assumption:

(θ| y) -(θ 0 | y) = 0,
where is the log-likelihood function, and thus by denition this is equivalent to

log f ∆X (1) (y 1 |θ) -log f ∆X (1) (y 1 |θ 0 ) + log f ∆X (2) (y 2 -(1 -ρ)y 1 |θ) f ∆X (2) (y 2 -(1 -ρ 0 )y 1 |θ 0 ) = 0. (2.11)
for all y 1 in (0, ∞) and y 2 in max (1 -ρ)y 1 , (1 -ρ 0 )y 1 , ∞ . This equation is assumed to be true, regardless of the observations set (y 1 , y 2 ). Considering that y 2 tends towards ∞, the quotient of the p.d.f.

tends towards 1 and we obtain that

log f ∆X (1) (y 1 |θ) -log f ∆X (1) (y 1 |θ 0 ) = 0.
Thus, Equation (2.11) entails that

log f ∆X (1) (y 1 |θ) -log f ∆X (1) (y 1 |θ 0 ) = 0 log f ∆X (2) (y 2 -(1 -ρ)y 1 |θ) -log f ∆X (2) (y 2 -(1 -ρ 0 )y 1 |θ 0 ) = 0 (2.12)
for all y 1 , y 2 , and where the rst equation depends on y 1 only while the second one depends on both y 1 and y 2 . Using the expression of f ∆X (1) , the rst equation in (2.12) becomes

log Γ α 0 T β0 b αT β Γ (αT β ) b α0T β 0 0 + αT β -α 0 T β0 log y 1 + (b 0 -b)y 1 = 0, (2.13)
which provides

         log Γ(α0T β 0 )b αT β Γ(αT β )b α 0 T β 0 0 = 0 αT β -α 0 T β0 = 0 b 0 -b = 0 (2.14)
Chapter 2. Parametric inference for the Arithmetic Reduction of Degradation models based on the fact that the three functions y 1 → 1, y 1 → y 1 and y 1 → log y 1 involved in (2.13) are linearly independent. Then, b = b 0 and αT β = α 0 T β0 . Now, in the same way, we can write the second equation in (2.12) as

log Γ α 0 T β0 (2 β0 -1) b αT β (2 β -1) Γ (αT β (2 β -1)) b α0T β 0 (2 β 0 -1) 0 + αT β (2 β -1) -1 log(y 2 -(1 -ρ)y 1 ) -α 0 T β0 (2 β0 -1) -1 log(y 2 -(1 -ρ 0 )y 1 ) + b 0 (y 2 -(1 -ρ 0 )y 1 ) -b(y 2 -(1 -ρ)y 1 ) = 0.
Then, replacing b by b 0 and αT β by α 0 T β0 , we obtain

log Γ αT β (2 β0 -1) b αT β (2 β -2 β 0 ) Γ (αT β (2 β -1)) + αT β (2 β -1) -1 log(y 2 -(1 -ρ)y 1 ) (2.15) -αT β (2 β0 -1) -1 log(y 2 -(1 -ρ 0 )y 1 ) + by 1 (ρ 0 -ρ) = 0.
We consider that y 1 tends towards 0, which leads to

log Γ αT β (2 β0 -1) b αT β (2 β -2 β 0 ) Γ (αT β (2 β -1)) + αT β (2 β -2 β0 ) log(y 2 ) = 0
and because the functions y 2 → 1 and y 2 → log y 2 are linearly independent, we have

αT β (2 β -2 β0 ) log(y 2 ) = 0
for all y 2 > 0, and thus β = β 0 . As a result, Equation (2.15) can be written as

αT β (2 β -1) -1 log y 2 -(1 -ρ)y 1 y 2 -(1 -ρ 0 )y 1 + by 1 (ρ 0 -ρ) = 0
which is equivalent to

αT β (2 β -1) -1 log 1 + ρ -ρ 0 y 2 /y 1 -(1 -ρ 0 ) + by 1 (ρ 0 -ρ) = 0.
Now, when y 2 tends towards ∞, we derive that by 1 (ρ 0 -ρ) = 0 for all y 1 > 0, hence ρ = ρ 0 . Merging this result with the previous, we have

           ρ = ρ 0 β = β 0 b = b 0 αT β = α 0 T β0
and nally θ = θ 0 .

The maximum likelihood estimate can be computed numerically by maximizing the log-likelihood function l (θ|y) = log L (θ|y) whose expression is

(θ|y) = s   α(nT ) β log b - n j=1 log Γ ∆A (j)   (2.16)
2.3. Maximum likelihood estimation

+ s i=1 n j=1 ∆A (j) -1 log g (i) j (ρ) -b s i=1 n j=1 g (i) j (ρ).
We now look for the zeros of the partial derivatives of l, and the equation

∂ b l (θ|y) = 0 implies that b = sα(nT ) β s i=1 n j=1 g (i) j (ρ)
. Now b is substituted by this expression into Equation (2.16), which provides the prole likelihood function l dened by

l ((α, β, ρ)|y) = s i=1 n j=1 ∆A (j) -1 log g (i) j (ρ) (2.17) -s       α(nT ) β       1 -log sα(nT ) β s i=1 n j=1 g (i) j (ρ)       - n j=1 log Γ ∆A (j)      
.

Finally, we obtain an estimate θ = (α, β, ρ, b) of θ given by:

θ =       arg max (α,β,ρ) l ((α, β, ρ)|y) , sα(nT ) β s i=1 n j=1 g (i) j (ρ)       . (2.18)

Application to the Arithmetic Reduction of Degradation model of order innity

Here the MLE method is studied in the case of a ARD ∞ model. However, the development of this method is exactly the same as in the previous section, the only dierence is the expression of Y jT with respect to ρ and the observations, and consequently the denition of the function g

(i) j .
For the ARD ∞ model, the degradation level Y jT right after the maintenance at time jT can be easily expressed with respect to Y j and the maintenance eciency parameter ρ, as we have

Y jT = (1 -ρ)Y j .
Then, we dene g

(i) j as g (i) j (ρ) = Y j -(1 -ρ)Y j-1 .
With such an expression of g (i) j , Lemma 1 holds for the ARD ∞ model and the likelihood function has the same expression provided by (2.16). As a result, the likelihood expression is the same here. Note that both models are identical over the time interval [0, 2T ), during which the rst two observations are conducted. Hence, the parameters of the ARD ∞ model are identiable as before. As a conclusion, Equations (2.17 For each parameters combination, 2000 data sets are generated, and an estimation of θ is computed for each data set. The ME method is here based on the minimization of a bivariate function, resulting in an estimation of β and ρ. The optimization is performed through a constrained gradient method: β is sought over the interval [0.1, 5] and ρ in [0.01, 0.99]. This iterative algorithm is initialized at point (2.5, 0.5) when ρ = 0.5, and at point (2.5, 0.8) otherwise. Regarding the MLE method, it is based on the maximization of a function depending on three parameters: α, β and ρ. Once again the optimization is performed through a constrained gradient method in the same way as before for β and ρ, and the parameter α is sought in the same interval as β. Some investigations lead us to conclude that the initialization point of the iterative optimization procedure has no inuence on the estimation results, and it only has very few impact on the CPU times. Moreover, the CPU time is negligible as the mean CPU time for one single estimation is less than 0.1 second for the ME method and 0.25 second for the MLE method.

The results are analysed through two indicators, the Relative Bias (RB) in percentage, dened as

RB θ = 100 × | θ -θ 0 | θ 0 ,
where θ is the empirical mean of the estimates, and the empirical variances of the estimates. For each possible combination of (β, ρ), the RB and the variances of the estimations are summarized in Appendix A. The gures deal either with the RB or with the variance, and are composed of four graphs: one for Chapter 3. Simulation study each parameter of the model. These graphs show four curves each, representing either the relative bias or the variance evolution with s for both the ME and the MLE methods when n = 5 and n = 20.

To be able to compare the results from one parameter to another, as well as the parameters sets, the y-axis range is the same for all plots. However, both the RB and the variance can be higher than the maximum value displayed in the y-axis. This signies that if a curve does not appear on a graph, as the black one in the top right-hand graph as well as the bottom left-hand one of Figure A.1, it means that all the points of the curve are out of range.

Analysis of the results

When n = 5, the ME method provides good estimations in terms of relative bias and variance of:

α and b when ρ = 0.2;

β when ρ = 0.5;

b when ρ = 0.8.

The value of β does not aect signicantly the estimations quality, but note that ρ is never well estimated here, even if the RB is sometimes low, the variance is always fairly high. As a result, this method provides poor quality estimations regarding the four parameters simultaneously.

If n = 20, once again the value of β does not aect signicantly the results. In this case, when ρ = 0.8, the ME method provides good estimations of the parameters. However, if ρ = 0.5 then the estimation quality of b is poor, it is also true for the estimation of ρ when its true value is set to 0.2.

Regarding the MLE method, it provides good quality estimations, better than those from the ME method. The RB of all the parameters estimators are less than 2%, except when ρ = 0.2, n = 5 and s = 50 where the RB related to ρ approaches 5%. Also, the variances are low, especially regarding β and ρ: they are between 5.4 × 10 -3 and 9.1 × 10 -5 for β, and between 2.2 × 10 -3 and 1.4 × 10 -6 for ρ. Here the value of β does not aect the results as well. Finally, the MLE method is more ecient in the case where n = 20 than when n = 5, except sometimes regarding the estimation of the scale parameter.

Remark 1. The case where β < 1 is not studied here due to numerical issues with the computation of the likelihood function given in Equation (2.10). Actually, if β < 1 then the shape parameter ∆A (j) of the degradation increments ∆X (j) is less than 1 as well. Hence the quantity g (i) j (ρ) ∆A (j) -1 might tend towards innity. To be more precise, there exists in some case a value ρ in (0, 1) such that the likelihood is only dened over (ρ, 1) regarding the parameter ρ, with g (i) j (ρ) = 0 and lim ρ→ ρ+ L θ| y (1) , . . . , y (s) = ∞.

Although this particularity allows to have a good estimation of the maintenance eciency (see Chapter 4), the sensitivity of the likelihood in the neighbourhood of ρ prevents the optimization from estimating properly the other parameters.

Arithmetic Reduction of Degradation model of order innity

Here the same analysis as in the previous section is done, with an identical framework for the model parameters, the observations characteristics, the optimization procedure as well as for the results summary which appears in Appendix B. Also, the case β < 1 is not studied for the reasons already mentioned 3.2. Arithmetic Reduction of Degradation model of order innity above.

Regarding the ME method, the parameter β is here well estimated when ρ = 0.2, associated with a low variance, and that whatever the values of n, s or β. When n = 5, this method is not reliable if ρ = 0.5 as none of the parameters is properly estimated, while α and b are rather well estimated if ρ = 0.8. When n = 20, the method provides good estimations of β and b if ρ = 0.5. Finally, as for the ARD 1 model, the ME method is ecient when ρ = 0.8, but the variance of this parameter is rather high. Once again the true value of β does not aect the results in a signicant way.

About the MLE method, the results are close to those of the previous section:

globally very good estimations with RB values smaller than 2%, except when ρ = 0.2;

when ρ = 0.2, the RB related to this parameter is higher (it approaches 3%);

the variances of the estimations of β and ρ are low, with similar boundaries as before;

this method is better as n increases, except sometimes for α.

As a result, for both ARD 1 and ARD ∞ models, the MLE method is the best. or nuclear power plants), which have lead to the development of the so-called reliability theory. For many years, only lifetime data were available and the rst reliability studies were focused on lifetime data analysis (see, e.g., [START_REF] Meeker | Statistical methods for reliability data[END_REF]), which still remains of interest in many cases. In that context and in case of repairable systems with instantaneous repairs, successive failure (or repair) times appear as the arrival points of a counting process, and failures hence correspond to recurrent events. As for the type of possible repairs, typical classical models are perfect (As-Good-As-New) and minimal (As-Bad-As-Old) repairs, leading to renewal and non homogeneous Poisson processes as underlying counting processes, respectively (see [START_REF] Barlow | Mathematical theory of reliability[END_REF]). The reality often lies in-between, leading to the class of imperfect repairs. Many models have been envisioned in the literature for their modeling, such as, e.g., virtual age models introduced by Kijima [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] and further studied in [START_REF] De Toledo | ARA and ARI imperfect repair models: Estimation, goodness-of-t and reliability prediction[END_REF]16], geometric processes [START_REF] Lam | The geometric process and its applications[END_REF] (extended in [START_REF] Bordes | Extended geometric processes: semiparametric estimation and application to reliability[END_REF]) or models based on reduction of failure intensity [START_REF] De Toledo | ARA and ARI imperfect repair models: Estimation, goodness-of-t and reliability prediction[END_REF]16]. See, e.g., [START_REF] Doyen | On geometric reduction of age or intensity models for imperfect maintenance[END_REF] for a recent account and extensions of such models. See also [START_REF] Pham | Imperfect maintenance[END_REF] for more references and other models.

Nowadays, the development of online monitoring and the increasing use of sensors for safety assessment make it possible to get specic information on the health of a system and on its eective evolution over time, without waiting for the system failure. This information is often synthesized into a scalar indicator, which can for instance stand for the length of a crack, the thickness of a cable, the intensity of vibrations, corrosion level, ... This scalar indicator can be considered as a measurement of the deterioration level Chapter 4. Semiparametric estimate of the maintenance actions eciency of the system. The evolution of this deterioration indicator over time is nowadays commonly modeled through a continuous-time and continuous-state stochastic process, which is often considered to have an increasing trend. Classical models include inverse Gaussian [START_REF] Ye | The Inverse Gaussian Process as a Degradation Model[END_REF] or Wiener processes with trend [START_REF] Hu | Optimum step-stress accelerated degradation test for wiener degradation process under constraints[END_REF][START_REF] Liu | A condition-based maintenance policy for degrading systems with age-and state-dependent operating cost[END_REF][START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF], which are also quite common in many other elds out of reliability theory, such as nance, insurance or epidemiology. This paper focuses on gamma processes, which are widely used since they were introduced in the reliability eld by Çinlar [START_REF] Çinlar | Stochastic process for extrapolating concrete creep[END_REF]. See [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF] and its references for a large overview.

In order to mitigate the degradation of the system over time and extends its lifetime, preventive maintenance actions can be considered, in addition to corrective repairs which are performed at failure.

In the context of deteriorating systems, many preventive maintenance policies from the literature consider condition-based maintenance (CBM) actions, where the preventive repair is triggered by the reaching of a preventive maintenance threshold by the deterioration level. In that context, "most of the existing CBM models have been limited to perfect maintenance actions", as noted by [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] (see also [START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF]). Some imperfect repair models are however emerging in the latest reliability literature, in this new context of deteriorating systems, see [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] for a recent review. Some models are based on the notion of virtual age previously introduced in the context of recurrent events (see, e.g., [START_REF] Giorgio | A new state-dependent degradation process and related model misidentication problems[END_REF][START_REF] Mercier | On the modelling of imperfect repairs for a continuously monitored gamma wear process through age reduction[END_REF]), where the system is rejuvenated by a maintenance action. Other models consider that an imperfect repair reduces the deterioration level of the system, such as [START_REF] Khatab | Condition-based selective maintenance for stochastically degrading multi-component systems under periodic inspection and imperfect maintenance[END_REF][START_REF] Letot | An adaptive degradation-based maintenance model taking into account both imperfect adjustments and agan replacements[END_REF][START_REF] Nicolai | Modelling and optimizing imperfect maintenance of coatings on steel structures[END_REF][START_REF] Ponchet | Maintenance policy on a nite time span for a gradually deteriorating system with imperfect improvements[END_REF][START_REF] Wu | A cost eective degradation-based maintenance strategy under imperfect repair[END_REF], which can be accompanied by some increase in the deterioration rate, as in [START_REF] Do | A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions[END_REF]. Also, some papers consider that the eciency decreases with the number of repair (see, e.g., [START_REF] Liu | An imperfect maintenance policy for mission-oriented systems subject to degradation and external shocks[END_REF][START_REF] Zhang | An ameliorated improvement factor model for imperfect maintenance and its goodness of t[END_REF]), and further studies, as in [START_REF] Huynh | Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems[END_REF], deal with imperfect maintenance models such that (i) repairs have a random eciency (ii) the deterioration rate increases with the number of repairs.

In all these papers however, the main point mostly is on the optimization of a maintenance policy, including these imperfect maintenance actions together with perfect repairs (replacements). Up to our knowledge, very few papers from the literature deal with statistical issues concerning imperfect repair models for deteriorating systems, except from [START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF], where the authors suggest a maximum likelihood method for estimating the parameters of the Wiener process together with an iterative procedure based on a Kalman lter for the dierent factors implied in successive imperfect repairs. This estimation procedure is developed in a fully parametric context and validated on simulated data, without any study of the asymptotic properties of the estimators.

The evaluation of the maintenance actions eciency is mainly used for maintenance policies optimization. Once the repair eciency has been estimated, the future behavior of the maintained system can be predicted, which allows to adapt (optimize) the periodicity of the maintenance actions and eciently plan a general overhaul. From a safety point of view, the principal inquiry is to ensure that the maintenance actions are eective enough to keep with a high probability the degradation level below a xed threshold (safety level). As long as this safety level is not reached, the maintenance actions must be adjusted, either by adapting their periodicity or by improving their eciency (if possible). Of course, apart from the previous safety concern, the maintenance costs are another issue. As an example, in [START_REF] Wu | A cost eective degradation-based maintenance strategy under imperfect repair[END_REF],

the costs minimization is based on the monitoring time and on the imperfect maintenance eciency. In [START_REF] Huynh | Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems[END_REF], the author considers a threshold for the degradation, beyond which an imperfect maintenance is performed. The optimization is made with respect to this threshold, the inspections periodicity and the repairs eciency. See, e.g., both papers cited above and their reference for an overview on maintenance policies optimization.

This paper focuses on a specic imperfect repair model, where each maintenance action reduces the deterioration level of the system. The model was rst introduced in [START_REF] Castro | Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs[END_REF] and further studied in [34], where it was called Arithmetic Reduction of Degradation model of order 1 (ARD 1 ). Mimicking Arithmetic Reduction of Age (ARA 1 ) and Arithmetic Reduction of Intensity (ARI 1 ) models of order 1 developed by 4.2. Framework [16] in the context of recurrent events, the idea of an ARD 1 model is that a maintenance action removes a proportion ρ of the degradation accumulated by the system from the last maintenance action (where ρ ∈ [0, 1)). The parameter ρ appears as a measure of the maintenance eciency, which lies between As-Good-As-New when ρ → 1 -and As-Bad-As-Old when ρ = 0. Along the same lines as [START_REF] Castro | Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs[END_REF]16,34], the maintenance actions eciency is here assumed to be xed and independent of the intrinsic degradation.

This paper is concerned with the development and study of an estimation procedure for the maintenance eciency parameter ρ, in the context of a gamma deteriorating system subject to periodic ARD 1 imperfect repairs. Observations are lead on just before each maintenance action. Considering n successive repairs, this leads to multivariate data, from where an estimator of ρ is proposed. The idea of this estimator has come from a preliminary study in a parametric framework based on the maximum likelihood method, where we have observed that the minimum of admissible ρ's has quite an interesting behavior, getting quickly very close to the unknown eciency parameter when n increases. This has lead to the proposition of an original estimator for ρ, which depends only on the data, and not on the shape function and rate parameter of the gamma process, leading to a semiparametric framework. Under technical assumptions, the strong consistency of this new estimator is shown, as the number n of repairs tends towards innity. Also, the convergence rate is proved to be surprisingly high, and can even reach an exponential speed in some cases. This estimator hence appears to be super consistent (under specic conditions). This is illustrated on simulated data at the end of the paper, where we provide two examples for which we observe that roughly 95% of the estimates are exact at the machine precision level (6.10 -17 )

as soon as n > 40 and n > 88, respectively, with a mean error below 10 -15 in both cases. The study is next extended to the case where s independent and identical systems are observed (n times each). A similar semiparametric estimator is proposed for the (common) maintenance eciency and the strong consistency is proved to hold as s tends towards innity, no matter the xed value of n and out of any technical condition requirement. The convergence rate is studied, which is shown to depend on the shape function of the gamma process and on the maintenance period, leading to a speed that can be either slower or faster than √ s, according to the case.

The outline of this paper is as follows. The framework is specied in Section 4.2, which covers the gamma deterioration process, the ARD 1 imperfect repair model and the observation scheme. Section 4.3

is devoted to the study of the semiparametric estimator in the case where one single system is observed, which includes its asymptotic properties when the number of repairs tends towards innity. Section 4.4 deals with the extension to several systems and considers the asymptotic properties with respect to the number of observed systems. Some illustrations of the estimator performances are provided in Section 4.5 and conclusions are formulated in Section 4.6.

Framework

Intrinsic deterioration

Let us rst recall that a random variable X is said to be gamma distributed with a and b as shape and rate parameters, respectively (X ∼ Γ(a, b) with a, b > 0), if its distribution admits the following probability density function (p.d.f.):

f X (x) = b a Γ(a) x a-1 e -bx 1 R + (x)
Chapter 4. Semiparametric estimate of the maintenance actions eciency with respect to Lebesgue measure. Its mean, variance and Laplace transform are provided by

E (X) = a b , V (X) = a b 2 , L X (t) = E e -tX = b b + t a , ∀t ≥ 0,
respectively. Moreover, cX ∼ Γ(a, b/c) for any c > 0, and the sum of n independent random variables

X i ∼ Γ(a i , b) (with 1 ≤ i ≤ n) is also gamma distributed with X 1 + • • • + X n ∼ Γ (a 1 + • • • + a n , b).
Now, let a(•) : R + → R + be a continuous and non decreasing function such that a(0) = 0 and let b > 0. Also let (X t ) t≥0 be a right-continuous process with left-side limits. Then, we recall that (X t ) t≥0 is a non homogeneous gamma process with shape function a(•) and rate parameter b, as soon as

X 0 = 0 almost surely (a.s.), (X t ) t≥0 has independent increments, each increment is gamma distributed: for all 0 ≤ s < t, we have X t -X s ∼ Γ(a(t) -a(s), b),
(see, e.g., [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF]).

In all the sequel, the intrinsic deterioration of the system (that is out of repairs) is assumed to be modeled by a non homogeneous gamma process (X t ) t≥0 with shape function a(•) and rate parameter b.

The imperfect repair model

In order to lower the deterioration level, instantaneous and periodic imperfect repairs are carried out on the system every T units of time (T > 0). Following [START_REF] Castro | Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs[END_REF]34], an Arithmetic Reduction Degradation model of order 1 (ARD 1 ) is considered, where a maintenance action removes a proportion ρ ∈ [0, 1) of the deterioration accumulated since the last maintenance action (or from time t = 0). The model used in the present paper is just the same as that used in [34], which we now recall, for sake of completeness.

Between repairs, the system is assumed to evolve according to independent and identically distributed (i.i.d.) copies X (i) t t≥0

, i = 1, 2, . . . of the gamma process (X t ) t≥0 , where exponent (i) refers to the i-th between-repair period [(i -1)T, iT ) (where time 0 is considered as a repair time). We set (Y t ) t≥0 to describe the overall deterioration level of the maintained system, as a result of the intrinsic deterioration and of the imperfect periodic repairs.

On the rst time interval [0, T ), there is no repair and

Y t = X (1) t for all t ∈ [0, T ). This implies that Y T -= X (1)
T -= X (1)
T a.s., based on the almost sure continuity of a gamma process. At time T , the deterioration level is reduced of ρX

(1) T so that Y T = (1 -ρ)X (1)
T a.s. On the second time interval [T, 2T ), we now have:

Y t = Y T + X (2) t -X (2)
T for all t ∈ [T, 2T ), which leads to

Y 2T -= Y T + X (2) 2T --X (2) T = (1 -ρ)X (1) T + X (2) 2T -X (2) T a.s.
and

Y 2T = (1 -ρ)X (1) T + (1 -ρ) X (2) 2T -X (2) T a.s.

Framework

More generally, on the j-th time interval (with j ∈ N), the eective degradation level can be expressed as

Y t = Y jT + X (j+1) t -X (j+1) jT for all t ∈ [jT, (j + 1)T ),
which leads to

Y jT -= (1 -ρ) j-1 p=1 X (p) pT -X (p) (p-1)T + X (j) jT -X (j) (j-1)T (4.1)
and

Y jT = (1 -ρ) j p=1 X (p) pT -X (p) (p-1)T
(with the convention that an empty sum is zero).

An example of trajectory of (Y t ) t≥0 is given in Note that out of maintenance times (t / ∈ {jT, j = 1, 2, . . . }), the random variable Y t is the sum of two gamma random variables which do not share the same rate parameter (except if ρ = 0). Hence, it is not gamma distributed. Please see [34] for more details on this model.

The periodicity T is assumed to be known and the previous model is called ARD 1 model with parameter (a(•), b, ρ) in the following (with T omitted).

As known from the introduction, our focus is on the development of an estimation procedure for the maintenance eciency parameter ρ. We now specify the observation scheme and derive some rst consequences.

Observation scheme and rst consequences

The deterioration level of the maintained system is assumed to be (perfectly) measured n times (n ∈ N * ), right before the n rst maintenance actions, that is at times T -, . . . , nT -. The data hence is an Chapter 4. Semiparametric estimate of the maintenance actions eciency observation of (Y T -, . . . , Y nT -), where Y jT -is provided by (4.1).

For the sake of simplicity, we set

Y j = Y jT -, U j = X (j) jT -X (j) (j-1)T, a j = a(jT ) -a (j -1)T for all j = 1, . . . , n and Y = (Y 1 , . . . , Y n ).
With the previous notation and based on the independent increments of a gamma process, the random variables U j 's can be seen to be independent with U j ∼ Γ(a j , b) for all j = 1, . . . , n. Also, for an ARD 1 model with parameter (a(•), b, ρ), Equation (4.1) can now be written as:

Y j = (1 -ρ) j-1 p=1 U p + U j (4.2)
for all 1 ≤ j ≤ n.

We rst check the theoretical identiability of the model, considering the parameters of the underlying gamma process as nuisance parameters (and T xed). Based on the independence between U 1 and U 2 , the Laplace transform of Y 2 is

L Y2 (t) = L (1-ρ)U1 (t) L U2 (t) = b b + (1 -ρ)t a1 b b + t a2 ,
with a similar expression for Ỹ2 . Remembering that a 1 = ã1 and b = b, this leads to

b b + (1 -ρ)t a1 b b + t a2 = b b + (1 -ρ)t a1 b b + t ã2 ,
for all t ≥ 0, which can be simplied into

1 + (1 -ρ)u 1 + (1 -ρ)u a1 = 1 (1 + u) ã2-a2
for all u ≥ 0, setting u = t/b. A rst order series expansion at point 0 induces that a 1 (ρ -ρ) = ã2 -a 2 and next that

1 + (1 -ρ)u 1 + (1 -ρ)u = 1 (1 + u) ρ-ρ 4.2.
Framework for all u ≥ 0. Taking the limit when u → +∞ in the previous equation, we get

1 - ρ 1 -ρ = lim u→+∞ 1 (1 + u) ρ-ρ =      0 if ρ > ρ, 1 if ρ = ρ, ∞ if ρ < ρ.
This is possible only if ρ = ρ since ρ and ρ belong to [0, 1), which achieves the proof.

The identiability hence holds as soon as two observations are available.

From now on, we assume that the true maintenance eciency parameter is ρ 0 ∈ [0, 1). The Y j 's and the U j 's hence correspond to an ARD 1 model with parameters (a (•) , b, ρ 0 ). A rst link between the Y j 's and the U j 's (j = 1, • • • , n) has been provided in Equation (4.2). We now invert this system of equations, thus providing an expression of the U j 's with respect to the Y j 's, that will be used in the sequel.

Lemma 2. For each j ∈ {1, • • • , n}, the increment U j can be expressed with respect to the observations Y 1 , . . . , Y j and to the maintenance eciency parameter ρ 0 as follows:

U j = j p=1 ρ j-p 0 (Y p -Y p-1 ) , (4.3)
where we set Y 0 = 0.

Proof. This result is proved by induction on j. For j = 1, the ARD 1 model denition provides

Y 1 = U 1 . Now, assume that (4.3) is true for some xed 1 ≤ j ≤ n -1. Observe from (4.2) that Y j+1 -Y j = (1 -ρ 0 ) j p=1 U p + U j+1 -(1 -ρ 0 ) j-1 p=1 U p + U j = U j+1 -ρ 0 U j , or equivalently that U j+1 = Y j+1 -Y j + ρ 0 U j .
Using the induction assumption, we easily derive that

U j+1 = Y j+1 -Y j + ρ 0 j p=1 ρ j-p 0 (Y p -Y p-1 ) = j+1 p=1 ρ j+1-p 0 (Y p -Y p-1 ) .
Hence, Equation (4.3) holds for j + 1, which achieves the proof.

For each 1 ≤ j ≤ n, let us now dene the function g j (ρ, Y) by

g j (ρ, Y) = j p=1 ρ j-p (Y p -Y p-1 ) , ∀ρ ∈ [0, 1), (4.4)
where we recall that the Y j 's refer to the true maintenance eciency parameter ρ 0 .

Lemma 2 ensures that g j (ρ 0 , Y) matches with the increment U j , that is for each 1 ≤ j ≤ n. As these increments are gamma distributed, they necessarily are non negative.

g j (ρ 0 , Y) = U j , (4.5
Hence the true parameter ρ 0 fulls the condition g j (ρ 0 , Y) ≥ 0 for each j ∈ {1, . . . , n}. An important consequence is that the range for the admissible ρ's can be restricted to the set Example 2. Two sets of parameters are considered, with ρ 0 = 0.5, T = 1 and b = 1 for both, a(t) = √ t (concave function) for the rst set and a(t) = t 1.5 (convex function) for the second one. An observation of Y is generated for each of the two parameter sets, and the corresponding observations of the D j 's are next computed. They are plotted in the left and right plots of Figure 4.2 for the concave and convex cases, respectively. The range for n is {1, • • • , 30} for the left plot (concave case) and 1, • • • , 10 6 for the right plot (convex case). Also, the parameter ρ 0 is highlighted by a vertical blue line on each plot. We can observe that in both cases, the sets D j 's are intervals of the shape [M j , 1) and that

D n = {ρ ∈ [0, 1) : g j (ρ, Y) ≥ 0 for all j ∈ {1, . . . , n}} .
M 1 ≤ M 2 ≤ • • • ≤ M n ≤ ρ 0 .
(Please note that the M j 's are indicated by blue crosses on the graphs). As can be seen on the left plot, it seems that, in case of a concave shape function, the sequence (M j ) converges very quickly towards ρ 0 when j increases. When the shape function is convex, it might also be convergent towards ρ 0 (?), but if so, it can only be at a very slow rate.

From the previous example, it seems that M n could be a very good estimator for ρ 0 in the concave case. However, in the convex case, even if the sequence (M j ) happened to converge towards ρ 0 when j increases (which we do not know), the rate of convergence would apparently be far below the classical square-root speed that could be obtained, e.g., with a maximum likelihood estimator. There hence is no interest in pursuing on this way in the convex case.

The semiparametric estimator and its asymptotic properties

As a summary, from the previous observations, we suggest to use M n as an estimator of the maintenance eciency parameter ρ 0 , that we hope to be convergent at a very high speed in the case of a concave shape function. Note that it is a semiparametric estimator of ρ 0 since the parameters of the gamma process are unknown and not restricted to a parametric family (the shape function a(•) is unknown).

The semiparametric estimator and its asymptotic properties

This section is devoted to the formal denition of the semiparametric estimator (Subsection 4.3.1), together with the study of its asymptotic properties (Subsection 4.3.3), when the number of imperfect repairs n tends to innity, with one single system observed.

Denition and rst properties

Let us rst recall from the previous section that

D j = {ρ ∈ [0, 1) : g k (ρ, Y) ≥ 0 for all k ∈ {1, . . . , j}} for all 1 ≤ j ≤ n. Note that each set D j is non empty, because g k (1 -, Y) = Y k ≥ 0 for all k ∈ {1, . . . , j}. Also, g 1 (ρ, Y) = Y 1 ≥ 0 for all ρ ∈ [0, 1), which implies that D 1 = [0, 1). Finally, it is readily seen that D j+1 = {ρ ∈ D j : g j+1 (ρ, Y) ≥ 0} and that D j+1 ⊂ D j for all 1 ≤ j ≤ n -1.
Let us set

M j = inf (D j )
for all 1 ≤ j ≤ n.

Proposition 7. The function g j (ρ, Y) and the sequence (M j ) 1≤j≤n (almost surely) satisfy the following properties:

1. M j ≤ M j+1 for all 1 ≤ j ≤ n -1; 2. ρ → g j (ρ, Y) is non negative on D j for all 1 ≤ j ≤ n; 3. ρ → g j (ρ, Y) is non decreasing on D j-1 for all 1 ≤ j ≤ n (where we set D 0 = D 1 = [0, 1) for convenience); 4. D j = [M j , 1) for all 1 ≤ j ≤ n; 5. M j ≤ ρ 0 for all 1 ≤ j ≤ n.
Proof. Points 1 and 2 are clear due to D j+1 ⊂ D j and to the denition of D j .

Let us show the three following points (Points 3-5) all together by induction on j.

At rst, we have

D 0 = D 1 = [0, 1), M 1 = 0 ≤ ρ 0 and g 1 (ρ, Y) = Y 1 for all ρ ∈ [0, 1
g j+1 (ρ, Y) = Y j+1 -Y j + ρ j p=1 ρ j-p (Y p -Y p-1 ) = Y j+1 -Y j + ρ g j (ρ, Y) (4.6)
for all 1 ≤ j ≤ n -1 (where Y j+1 -Y j might be negative).

By the induction assumption, g j (ρ, Y) is non decreasing on D j-1 , and hence also on D j (as D j ⊂ D j-1 ). Based on the previous recursion formula (4.6) and as g j (ρ, Y) also is non negative on D j , this implies that g j+1 (ρ, Y) is non decreasing on D j .

As D j = [M j , 1) by the induction assumption, we now have

M j+1 = inf {ρ ∈ [M j , 1) : g j+1 (ρ, Y) ≥ 0} ,
where g j+1 is non decreasing and continuous on [M j , 1),

with g j+1 (1 -, Y) = Y j+1 > 0 (almost surely).
This implies that D j+1 is an interval and D j+1 = [M j+1 , 1).

Finally, from Equation (4.5), we have g j+1 (ρ 0 , Y) = U j+1 ≥ 0. As M j ≤ ρ 0 by the induction assumption and g j+1 (ρ, Y) is known to be non decreasing on [M j , 1), we necessarily have M j+1 ≤ ρ 0 .

Hence Points 2-5 are true for j + 1, and this achieves the proof.

Based on the previous results, we can see that the sequence (M j ) 1≤j≤n can be alternatively dened through

M 1 = 0, M j+1 = inf {ρ ∈ [M j , 1) : g j+1 (ρ, Y) ≥ 0} , for all 1 ≤ j ≤ n -1. ( 4

Technical results

Lemma 3. Let ρ ∈ [0, 1). Then,

g j (ρ, Y) ≥ 0 implies that ρ 0 -ρ ≤ U j U j-1 for each 2 ≤ j ≤ n.
Proof. Let us rst prove by induction that

g j (ρ, Y) = (ρ -ρ 0 ) j-1 p=1 ρ j-1-p U p + U j • For j = 1, the result is clear because g 1 (ρ, Y) = Y 1 = U 1 .
Assume that the result holds for some

j ∈ {1, • • • , n -1}.
Based on (4.6) and (4.5), we know that

g j+1 (ρ, Y) = Y j+1 -Y j + ρ g j (ρ, Y) , U j+1 = Y j+1 -Y j + ρ 0 U j ,
4.3. The semiparametric estimator and its asymptotic properties (taking ρ = ρ 0 in the rst line to derive the second one), which provides

g j+1 (ρ, Y) = U j+1 -ρ 0 U j + ρ g j (ρ, Y) .
Using the induction assumption, g j+1 (ρ, Y) can now be expressed as follows:

g j+1 (ρ, Y) = U j+1 -ρ 0 U j + ρ (ρ -ρ 0 ) j-1 p=1 ρ j-1-p U p + U j = (ρ -ρ 0 ) j p=1 ρ j-p U p + U j+1
where the last equality results from straightforward calculations. Thus we obtain the rst result.

Next we note that g j (ρ, Y) ≥ 0 is true as soon as

(ρ -ρ 0 ) j-1 p=1 ρ j-1-p U p + U j ≥ 0,
or equivalently

ρ 0 -ρ ≤ U j U j-1 + j-2 p=1 ρ j-1-p U p .
This implies the result since

j-2 p=1 ρ j-1-p U p ≥ 0.
In the following, we will have to control quantities of the shape P (ρ 0 -ρ > ε), which will be done by controlling quantities of the shape P (U j /U j-1 > ε), using arguments based on the previous lemma. This will be achieved through the use of the following Remark and Lemma.

Remark 2. For each j ≥ 2, the random variables U j-1 and U j are known to be independent and gamma distributed Γ (a j-1 , b) and Γ (a j , b), respectively. It follows that, for all ε ≥ 0,

P U j U j-1 > ε = P U j U j-1 + U j > ε 1 + ε
where the random variable U j / (U j-1 + U j ) is beta distributed B (a j , a j-1 ) (standard property of gamma distributions), with p.d.f. with positive parameters α 1 and α 2 (which is also called the regularized incomplete beta function). For all x ∈ [0, 1],

f aj ,aj-1 (t) = 1 B(a j , a j-1 ) t aj -1 (1 -t) aj-1-1 , ∀t ∈ [0, 1] .
I x (α 1 , α 2 ) ≥ x α1 (1 -x) α2 1 + α1 α2 . (4.9)
Chapter 4. Semiparametric estimate of the maintenance actions eciency

Proof. Let us rst show that

I x (α 1 , α 2 ) ≥ x α1 (1 -x) α2 α 1 B(α 1 , α 2 ) (4.10) for all x ∈ [0, 1].
Note that Inequality (4.10) can be seen as a direct consequence of [14,Eq. 8.17.20]. As it is stated without proof in the quoted reference, we prefer to propose some details here.

For (α 1, α 2 ) xed, let us set

g (x) = I x (α 1 , α 2 ) - x α1 (1 -x) α2 α 1 B(α 1 , α 2 ) , ∀x ∈ [0, 1] .
Based on the p.d.f. of a beta distribution recalled in (4.8), it is easy to check that

g (x) = x α1 (1 -x) α2-1 B(α 1 , α 2 ) 1 + α 2 α 1 ≥ 0.
Thus g (x) is non decreasing with respect to x. As g (0) = 0, we derive that g (x) ≥ 0 for all x ∈ [0, 1]

and Inequality (4.10) is true.

It remains to show that

1 α 1 B(α 1 , α 2 ) ≥ 1 1 + α1 α2 (4.11)
to derive (4.9). Now, [2, Eq. 6.1.3 p. 255] states that for all positive real number α, the inverse of Γ(α) can be expressed as

Γ(α) -1 = α exp (γα) m≥1 1 + α m exp - α m ,
where γ is Euler's constant. By denition of the Beta function, we hence have

B(α 1 , α 2 ) = Γ (α 1 ) Γ (α 2 ) Γ (α 1 + α 2 ) = α 1 + α 2 α 1 α 2    m≥1 1 + α1+α2 m exp -α1+α2 m m≥1 1 + α1 m exp -α1 m m≥1 1 + α2 m exp -α2 m    .
As the products are convergent, this can be simplied into

B(α 1 , α 2 ) = α 1 + α 2 α 1 α 2 m≥1 1 + α1+α2 m 1 + α1 m 1 + α2 m ≤ 1 α 1 + 1 α 2 ,
which provides (4.11) and the result.

Corollary 1. Let 2 ≤ j ≤ n. Then ρ 0 -M j ≤ U j /U j-1 and

P (ρ 0 -M j > ε) ≤ P U j U j-1 > ε ≤ 1 - εaj (1 -ε) aj-1 1 + aj aj-1
4.3. The semiparametric estimator and its asymptotic properties for all ε ∈ (0, 1), with ε = ε/ (1 + ε).

Proof. By denition of M j , we know that g j (M j , Y) ≥ 0. Based on Lemma 3, we derive that ρ 0 -M j ≤ U j /U j-1 . Hence:

P (ρ 0 -M j > ε) ≤ P U j U j-1 > ε .
Now, a direct consequence from Remark 2 and Lemma 4 can be written as follows:

P U j U j-1 > ε = 1 -I ε(a j , a j-1 ) ≤ 1 - εaj (1 -ε) aj-1 1 + aj aj-1 with ε = ε/ (1 + ε), which allows to conclude.
All the previous results are valid without any assumption on the shape function a (•). We now come to specic technical results, which requires some concavity assumption for a (•) to hold. To be more precise, our main assumption states as follows.

Assumption (H) The shape function a (•) of the gamma process is concave and dierentiable on R + , and such that lim t→∞ a (t) = 0. Remark 3. All the asymptotic results of the paper remain valid if the concavity and dierentiability properties for a (•) only hold from one point t 0 (that is on a set [t 0 , +∞)), and not on the whole half-line R + . Remark 4. Many classical concave shape functions fullls Assumption (H), such as

a 1 (t) = αt β , a 2 (t) = log 1 + αt β or a 3 (t) = 1 -exp -αt β ,
with α > 0 and 0 < β < 1 for the rst case and 0 < β ≤ 1 in the other two cases.

When β > 1, both shape functions a 2 (•) and a 3 (•) are concave only from the point t 0 , with t 0 = [(β -1)/α] 1/β and t 0 = [(β -1)/(αβ)] 1/β , respectively. Hence, as stated in Remark 3, all the asymptotic results remain valid for a (•) = a 2 (•) or a 3 (•) and for all β > 0.

Nevertheless, for the sake of simplicity, the shape function a (•) is assumed to be concave and differentiable from the initial time, in all the results requiring Assumption (H) to hold. Lemma 5. Suppose Assumption (H) to hold. Then the sequence (a n ) n∈N * is non increasing, and tends to 0 when n tends to ∞.

Proof. By the mean value theorem, there exists c(n) in ((n -1)T , nT ) such that

a n = a(nT ) -a((n -1)T ) = T a (c(n)).
As lim t→∞ a (t) = 0 by assumption and lim n→∞ c(n) = ∞, it follows that lim t→∞ a (c(n)) = 0, which induces the convergence of (a n ) n∈N * towards 0. Finally, a n ≥ a n+1 for all n ≥ 0 is a direct consequence of the concavity of a(•). Lemma 6. Suppose Assumption (H) to hold. Then,

P (ρ 0 -M 4n > ε) ≤ 1 - for all n ≥ 1 and all ε ∈ (0, 1). Proof. Let 2 ≤ k ≤ n. Based on Corollary 1, we know that ρ 0 -M k ≤ U k /U k-1 . As (M k ) 2≤k≤n is non decreasing, we get that ρ 0 -M n = min 2≤k≤n (ρ 0 -M k ) ≤ min 2≤k≤n U k U k-1 .
Now, in order to boil down to independent random variables U k /U k-1 , let us consider only the even terms k = 2j. Also, for sake of simplication, let us substitute n by 4n. We get

ρ 0 -M 4n ≤ min 2≤k≤4n U k U k-1 ≤ min 1≤j≤2n U 2j U 2j-1 .
This induces

P (ρ 0 -M 4n > ε) ≤ P min 1≤j≤2n U 2j U 2j-1 > ε = 2n j=1 P U 2j U 2j-1 > ε
since the ratios are independent random variables. Based on Corollary 1 again, we derive that

P (ρ 0 -M 4n > ε) ≤ 2n j=1 1 - ε a2j (1 -ε) a2j-1 1 + a2j a2j-1 (4.12)
where ε = ε/ (1 + ε). Under Assumption (H), we know from Lemma 5 that a 2j ≤ a 2j-1 for all 1 ≤ j ≤ 2n.

Thus

1 1 + a2j a2j-1 ≥ 1 2 and 2n j=1 1 - ε a2j (1 -ε) a2j-1 1 + a2j a2j-1 ≤ 2n j=1 1 - 1 2 ε a2j (1 -ε) a2j-1 ≤ 2n j=n+1 1 - 1 2 ε a2j (1 -ε) a2j-1 (4.13)
(as each term is smaller than 1 in the product).

Using again the non increasingness of (a j ) n+1≤j≤2n , we get that εa2j ≥ εa 2(n+1) and

(1 -ε) a2j-1 ≥ (1 -ε) a2n+1 for all j ∈ {n + 1, n + 2, . . . , 2n}, since ε ∈ (0, 1). This implies 2n j=n+1 1 - 1 2 ε a2j (1 -ε) a2j-1 ≤ 1 - 1 2 ε a 2(n+1) (1 -ε) a2n+1 n . (4.14)
Putting together (4.12), (4.13) and (4.14) leads to

P (ρ 0 -M 4n > ε) ≤ 1 - 1 2 ε a 2(n+1) (1 -ε) a2n+1 n = 1 - ε a 2(n+1) 2 (1 + ε) a2n+1+a 2(n+1) n 4.3.
The semiparametric estimator and its asymptotic properties by denition of ε. Finally, because 1/ (1 + ε) > 1/2 (as ε < 1), we have

P (ρ 0 -M 4n > ε) ≤ 1 - ε a 2(n+1) 2 1+a2n+1+a 2(n+1)
n , which nishes the proof.

We are now ready to state our main results, which is done in next subsection.

Consistency and convergence rates

Theorem 1. Suppose Assumption (H) to hold. Then M n is a strongly consistent estimator of ρ 0 (M n -→ ρ 0 almost surely) as the number of repairs n tends to innity.

Proof. First, because (M n ) n≥1 is non decreasing, it is sucient to prove the almost sure convergence of the subsequence (M 4n ) n≥1 towards ρ 0 . (The same remark is valid for the convergence rate, hereafter).

Let ε ∈ (0, 1). From Lemma 6, we know that n≥1

P (ρ 0 -M 4n > ε) ≤ n≥1 1 - ε a 2(n+1) 2 1+a2n+1+a 2(n+1) n . (4.15)
Hence, it is enough to show the convergence of the right-side series in the previous inequality, to show the strong consistency.

Under Assumption (H) and by Lemma 5, we have lim

n→+∞ a 2(n+1) = lim n→+∞ a 2n+1 = 0. Then lim n→+∞ n 1 - ε a 2(n+1) 2 1+a2n+1+a 2(n+1) n = lim n→+∞ 1 - ε a 2(n+1) 2 1+a2n+1+a 2(n+1) = 1 2 < 1.
The root test ensures the convergence of the right-side series in (4.15), which allows to conclude.

We now look at the convergence rate, which reveals itself to be very high (at least sub-exponential, or even exponential).

Theorem 2. Suppose Assumption (H) to hold. Then we have:

1. The almost sure convergence rate of the estimator M n is at least sub-exponential (that is at least polynomial of order k, for any k > 0) as soon as

a 2n = O (log n) -1 .
2. The almost sure convergence rate is at least exponential as soon as

a 2n = O n -1 .
3. The convergence rate in probability is at least exponential as soon as a 2n = o n -1 log n .

Proof. Let ε n ∈ (0, 1) for all n ∈ N * . Based on Lemma 6, we have

P (ρ 0 -M 4n > ε 4n ) ≤ 1 - ε a 2(n+1) 4n 2 1+a2n+1+a 2(n+1) n .
As a 2n+1 ≤ a 2(n+1) ≤ a 1 from Lemma 5, we get that 2 1+a2n+1+a 2(n+1) ≤ 2 1+2a1 and hence and C = 1/ 2 1+2a1 . From the root test, we know that the series with generic term (1 -u n ) n is convergent as soon as lim sup n→+∞ (1 -u n ) < 1, or equivalently as soon as lim inf n→∞ u n > 0. Hence, the series with generic term the left-side expression in (4. 16) is convergent as soon as lim inf n→∞ u n > 0.

P (ρ 0 -M 4n > ε 4n ) ≤ (1 -u n ) n
Let us now look at the three dierent points of the theorem. Assume rst that ε n = ε/n k with k > 0 and ε ∈ (0, 1). This provides

u n = C ε (4n) k a 2(n+1) = C exp a 2(n+1) [log (ε) -k log (4) -k log (n)] .
Assume further that a 2n = O (log n) -1 , or equivalently that a 2(n+1) = O (log n) -1 . Then, there exists K > 0 such that a 2(n+1) log (n) < K, from where we derive that

u n > C exp a 2(n+1) (log (ε) -k log (4)) -k K .
Hence

lim inf n→∞ u n ≥ C lim inf n→∞ exp {-k K} > 0
because a 2(n+1) converges towards 0 (see Lemma 5).

This shows that the series with generic term the left-side expression in (4. 16) is convergent for ε n = ε/n k and any (k, ε), which means that M n almost surely converges towards ρ 0 at speed at least n -k for any k > 0, namely the convergence rate is at least sub-exponential, which proves the rst point. Now let us set ε n = ε exp(-kn) with k > 0 and ε ∈ (0, 1).

We have

u n = C ε a 2(n+1) exp -4kna 2(n+1) = C exp na 2(n+1) log (ε) n -4k .
Assume that a 2n = O n -1 . Then, there exists K > 0 such that na 2(n+1) < K, from where we derive that

u n > C exp K log (ε) n -4k . Hence lim inf n→∞ u n ≥ C exp (-4Kk) > 0,
which allows to conclude for the second point.

Finally, assume that a 2n = o n -1 log n . The point here is to show the convergence in probability.

Based on (4.16), it is sucient to show that lim n→+∞ (1 -u n ) n = 0.

We have

(1 -u n ) n = exp n log 1 -C ε a 2(n+1) exp -4kna 2(n+1) .

The semiparametric estimator and its asymptotic properties

As log (1 -x) ≤ -x for all x ∈ (0, 1), we get that

(1 -u n ) n ≤ exp -C nε a 2(n+1) exp -4kna 2(n+1) = exp (-C v n ) (4.17) with v n = exp log (n) 1 + na 2(n+1) log (n) log (ε) n -4k
Based on a 2n = o n -1 log n , we have lim n→+∞ na 2n / log (n) = 0, which implies

lim n→+∞ 1 + na 2(n+1) log (n) log (ε) n -4k = 1
and hence v n tends to ∞. We derive from (4.17) that (1 -u n ) n converges towards 0, which allows to conclude.

Example 3. Let a(t) = αt β with 0 < α, β < 1, which is already known to fulll Assumption (H) from Remark 4. Also, we have

a 2n = αT β (2n) β -(2n -1) β = α(2nT ) β 1 -1 - 1 2n β ∼ n→+∞ C 1 n 1-β
where C = αβ2 β-1 T β . It is easy to check that the condition a 2n = O (log n) -1 from Point 1 in Theorem 2 is satised (but not the conditions for the other points). Hence, we can conclude that the almost sure convergence holds with an at least sub-exponential rate.

Example 4. Let a(t) = log 1 + αt β with α > 0, 0 < β ≤ 1, which is already known to fulll Assumption (H) from Remark 4. Also, based on log (x) ∼ x -1when x → 1 for the second line, we have

a 2n = log 1 + αT β (2n) β 1 + αT β (2n -1) β ∼ n→+∞ αT β (2n) β -(2n -1) β 1 + αT β (2n -1) β ∼ n→+∞ 1 (1 -1 2n ) β -1 ∼ n→+∞ β 2n .
Hence, the condition a 2n = O n -1 is satised (strongest condition in Theorem 2), and the almost sure convergence holds with an at least exponential rate. Note that this result would remain valid for β > 1 as the shape function is concave from point t 0 = [(β -1)/α] 1/β (see Remark 4).

Example 5. Let a(t) = 1 -exp -αt β with α > 0, 0 < β ≤ 1, which is already known to fulll Assumption (H) from Remark 4. We have

a 2n = exp -α (2nT ) β -exp -α (2n -1) T β ,
which clearly implies a 2n = O n -1 . Hence the almost sure convergence holds with an at least exponential rate. Note that, here again, the result would remain valid for β > 1 as the shape function is concave from point t 0 = [(β -1)/(αβ)] 1/β (see Remark 4).

Up to here, it was supposed that one single system is observed. In the next section, we now envision the possibility of observing several systems. 2T -, ..., nT -, as described in Section 4.2. For each i ∈ {1, . . . , s}, we add exponent (i) to each quantity referring to the i-th system. For instance,

Y (i) = Y (i) 1 , . . . , Y (i) n
stands for the multivariate observation of the i-th system at times T -, 2T -, ..., nT -. Also, the sequence

M (i) j 1≤j≤n
is dened by

M (i) 1 = 0 M (i) j = inf ρ ∈ M (i) j-1 , 1 : g j ρ, Y (i) ≥ 0 for all 2 ≤ j ≤ n
in a similar way as in (4.7).

The extended semiparametric estimator is dened as

M s,n = max 1≤i≤s M (i)
n for all n ≤ 1 and s ≥ 1.

The asymptotic properties of each sequence

M (i) n n∈N *
(with i xed) has been studied in the previous section. Clearly, similar results are valid for the sequence (M s,n ) n∈N * with s xed (with an even higher rate of convergence as

M (i) n ≤ M s,n ≤ ρ 0 for each i).
We hence focus on the asymptotic properties of (M s,n ) s∈N * with n xed in the sequel of this section. We take n ≥ 2, which ensures the identiability, based on Proposition 6.

Consistency and convergence rates according to the number of observed systems

Theorem 3. Let n ≥ 2. Then M s,n is a strongly consistent estimator of ρ 0 as the number of observed systems s tends to innity.

Proof. From Proposition 7, we know that

M (i) 2 ≤ M (i) n ≤ ρ 0 for each 1 ≤ i ≤ s, which implies that M s,2 ≤ M s,n ≤ ρ 0 .
Hence, it is enough to prove that M s,2 is a strongly consistent estimator of ρ 0 . From the denition of M s,2 , we have

ρ 0 -M s,2 = min 1≤i≤s ρ 0 -M (i) 2
and as the systems are i.i.d., this provides

P (ρ 0 -M s,2 > ε) = P ρ 0 -M (1) 2 > ε s (4.18)
4.4. Extension to the case where several systems are observed for all ε > 0. Now, Corollary 1 leads to Proof. Let us set ε s = ε/s k , with k > 0 and ε ∈ (0, 1). The point is to show that the series with generic term P (ρ 0 -M s,n > ε s ) converges for all k < 1/ min 2≤p≤n a p . Using a similar procedure as for Equation (4.18) and based on Corollary 1, we have 2

P (ρ 0 -M s,2 > ε) ≤ P U (1) 2 U (1) 1 > ε s and because P U (1) 2 /U (1) 1 > ε < 1, we get s≥1 P U (1) 2 U (1) 1 > ε s < ∞.
P (ρ 0 -M s,n > ε s ) = P min 1≤i≤s ρ 0 -M (i) n > ε s ≤ n p=2 P p (s) s (4.19)
where

P p (s) = 1 - ε ap s 1 + ap ap-1 (1 + ε s ) ap-1+ap
.

Note that P p (s) ∈ (0, 1) for any s ≥ 1 and 2 ≤ p ≤ n so that the product in (4. [START_REF] Giorgio | A new state-dependent degradation process and related model misidentication problems[END_REF]) is smaller than each of its term. Hence, keeping only the 2p-th term, we get

P (ρ 0 -M s,2n > ε s ) ≤ (P 2p (s)) s
for all p ≥ 1 such that 2p ≤ n. Then, using that ε s = εs -k and 1/(1 + ε s ) > 1/2, we obtain

(P 2p (s)) s ≤   1 - s -ka2p 1 + a2p a2p-1 2 a2p-1+a2p ε a2p   s .
Now, using that log(1 -x) ≤ -x for all x in [0, 1), it follows that (P 2p (s)) s ≤ u s , where

u s = exp -C p s 1-ka2p and C p = ε a2p 1 + a 2p a 2p-1 2 a2p-1+a2p -1
.

Gathering the previous inequalities, we now have P (ρ 0 -M s,n > ε s ) ≤ u s and the point is to study the 2 To be more precise, we have P min s 2 , which entails that the series with generic term u s is convergent, and hence also the series with generic term P (ρ 0 -M s,n > ε s ). This allows to conclude that the almost sure convergence rate of the estimator M s,n is at least s -k , for any k < a -1 2p and any p ≥ 1 such that 2p ≤ n. Now, keeping only the 2p + 1-th term in the product in (4. [START_REF] Giorgio | A new state-dependent degradation process and related model misidentication problems[END_REF] provides

P (ρ 0 -M s,2n > ε s ) ≤ (P 2p+1 (s)) s
for any p such that 2p + 1 ≤ n. Similar arguments as above allow to derive that the convergence rate it at least s -k , for any k < a -1 2p+1 and any p ≥ 1 such that 2p + 1 ≤ n. Finally, it is hence true for any k < a -1 p and any p such that 2 ≤ p ≤ n, and consequently for any k < max 2≤p≤n a -1 p , which allows to conclude.

The real number a p corresponds to the increment of the shape function over the time interval [(p -1)T, pT ).

The overall convergence rate obtained in the previous theorem corresponds to the smallest increment.

Hence the smaller this increment, the slower the degradation and the higher the convergence rate. More precisely, when the shape function is concave, the increments decrease over time and the convergence rate is at least s -1/an because the smallest increment is the last one. On the other hand, the increments increase over time when the shape function is convex, hence the smallest increment is the second one and the convergence rate is at least s -1/a2 . Then the convergence rate is higher than the standard square-root speed as soon as the smallest increment is less than 2. Note that this condition depends on both the shape function and the period T , as illustrated in the next example.

Example 6. Let a(t) = αt with α > 0, hence a p = αT for each 2 ≤ p ≤ n and the almost sure convergence rate of M s,n with respect to s is at least s -1/αT . In comparison with the classical rate s -1/2 , it is higher if αT < 2 and lower if αT > 2.

Empirical illustration based on simulated data

The aim of this section is to illustrate our most signicant results, that is the fast convergence rates obtained in Section 4.3 in the case where a single system is observed. In that case asymptotic results are obtained with respect to an increasing number of repairs n. The point hence is to observe the empirical behavior of the semiparametric estimator M n of ρ 0 , which from Theorems 1 and 2 is known to be strongly consistent, with a convergence rate that can be either exponential (ECR) or sub-exponential (S-ECR)

with respect to n (considering either almost sure convergence or convergence in probability).

These illustrations are based on two simulated datasets, considering a(t) = log(1 + t) (rst case) and a(t) = √ t (second case) as shape functions, which from Examples 4 and 3 provide exponential and sub-exponential almost sure convergence rates respectively. The rst (resp. second) case will hence be referred to as ECR (resp. S-ECR) in the sequel.

Empirical illustration based on simulated data

Data simulation and empirical bias

To be able to compare results, we place ourselves within the same framework for both cases. First, the model parameters as well as the observation characteristics of a maintained system are the following: Chapter 4. Semiparametric estimate of the maintenance actions eciency Our point now is to explain these two particularities, which are induced by computing limitations, as is discussed hereafter.

Let M num be the largest positive real number such that 1 -M num = 1 (numerically negligible), which is roughly equal to 1.11 × 10 -16 in our case. This entails that (1 -M num )/2 = 1/2, and hence

ρ 0 - M num 2 = ρ 0 ,
based on ρ 0 = 1/2. Then, when we obtain an estimate m (i) n = 0.5 for some 1 ≤ i ≤ 250 000 and 1 ≤ n ≤ 250, it only means that ρ 0 -m (i) n ≤ M num /2 and not that m (i) n = ρ 0 . Thus, the numerical estimate of the bias is correct only if ρ 0 -m (i) n ≥ M num /2, otherwise it is underestimated. This leads us to introduce the Proportion of Numerically Exact Estimates (PNEE) as a function of n As a summary, both jumps and convergence to a non-zero constant in Figure 4.3 are due to numerical limitations and the corresponding points on the plots should not be considered for the further study of the empirical bias. Thus, in the sequel, we only focus on the rst values of n for which there is no jump, that is on 1 ≤ n ≤ 26 for the ECR case and 1 ≤ n ≤ 49 for the S-ECR case.

PNEE(n) = # 1 ≤ i ≤ 250 000 : ρ 0 -m (i) n = Mnum
Our aim now is to explore the convergence rate from an empirical point of view. This can be done through the study of the empirical bias, as is now explained.

Link between the bias and the exponential convergence rate

Let us show that the convergence rate of M n is at least exponential whenever log 10 ρ 0 -E (M n ) decreases linearly. Indeed, assuming that log 10 ρ 0 -E (M n ) decreases linearly, there exist k > 0 and C ∈ R such that

log 10 ρ 0 -E (M n ) = C -kn.
Because M n ∈ [0, ρ 0 ] for any n ≥ 1 with probability one, the random variable ρ 0 -M n is non negative and by the Markov's inequality we have

P ρ 0 -M n > ε n ≤ ρ 0 -E (M n ) ε n • Setting ε n = ε exp(-kn) with k ∈ (0, k) and ε > 0, we have P ρ 0 -M n > ε n ≤ exp(C) ε exp -k -k n -→ n→∞ 0,
and n≥1 P ρ 0 -M n > ε n ≤ n≥1 exp(C) ε exp -k -k n < +∞.
Therefore the two last results allow to conclude that the rate of convergence of M n is at least exponential for the convergence in probability as well as for the almost sure convergence. If

log 10 ρ 0 -E (M n )
decreases at a slower rate than the linear rate, we can not conclude to a sub-exponential convergence rate, however we observe that the empirical evidence of a slower convergence rate for log 10 ρ 0 -E (M n )

coincides with a slower convergence rate for ρ 0 -M n in Theorem 2.

Link between theoretical and empirical results

Because of the previous explanations on the behavior of EB(n), linear regressions are performed on {EB(n); 1 ≤ n ≤ 26} for the ECR case and {EB(n); 1 ≤ n ≤ 49} for the S-ECR case. We assume that the relationship between EB(n) and n is modelled by a linear regression either simple or quadratic. The results are summarized on Figure 4.5, Table 4.1 and Table 4.2.

The tables provide the linear regressions outcomes, that is the coecients of the rst and second degree polynomial regression (and their related p-values and indicators about residual error (minimum, maximum, rst, second and third quartiles). For the ECR case, both rst (simple) and second (quadratic) degree polynomial regression t well log 10 (EB(n)), with adjusted R 2 of 0.9994 and 0.9996, respectively.

However, Table 4.1 shows that the second degree term is not signicant neither in comparison with Regarding the optimality of our results, we recall that the condition of Theorem 2 to obtain an at least exponential convergence rate is a 2n = o n -1 log n . Now repeating the study with a 2n = n -1 log n we see in Figure 4.6 that again an exponential convergence rate is expected because the bias decreases linearly.

We conclude that the condition on a 2n in Theorem 2 is probably sucient but not necessary.

Concluding remarks

In this paper we propose a semiparametric inference method for the maintenance eciency parameter involved in the ARD 1 repair model for a Gamma deteriorating system. For a single system the main condition that insures the strong consistency of our semiparametric estimator of the maintenance eciency parameter is the concavity of the shape function of the underlying Gamma deteriorating process. Two types of asymptotic results are obtained: either a single system is observed with the number of repairs tending to innity, or it is the number of systems that tends to innity. In the case of a single system the almost sure convergence rate of the estimator can be particularly fast, at least exponential for some particular cases. The simulation study illustrates the convergence rates obtained in case of a single trajectory. We observe that the theoretical convergence rates are consistent with the numerical simulation results. However it seems that it is still possible to rene the mathematical conditions under which the convergence rates are obtained. Thus improving the assumptions accuracy may constitute further work.

Note that when several systems are considered the convergence rate of the estimator is slower but its strong consistency holds whatever the shape function is. Depending on the shape function, the convergence rate of the estimator may overcome the usual square root rate.

In the wake of this study we want to mention several lines of work that we consider as are important.

First, the observation scheme could be decoupled from the scheme of maintenance actions. Indeed, it would be interesting, for instance for an ARD 1 model, to consider a system with scheduled times of maintenance actions for which the observation times are independent of the maintenance schedule. Also, Chapter 4. Semiparametric estimate of the maintenance actions eciency there exist many other models that extend the ARD 1 model, such as for instance the ARD m (resp.

ARD ∞ ) for which the basic idea is that a maintenance action removes a proportion of the degradation accumulated by the system from the last m maintenance actions (resp. since the system was put into operation).

As explained in the introduction, there exist also alternatives to arithmetic reduction of degradation models such as those based on arithmetic reduction of age. For such models, instead of reducing the degradation level of the system, the maintenance action consists in reducing the age of the system. The use of these models is not restricted to gamma processes, and can be generalized to any non homogeneous Lévy process. As an example, [START_REF] Kahle | Imperfect repair in degradation processes: A kijima-type approach[END_REF] deals with both the ARD 1 and ARD ∞ models, as well as two arithmetic reduction of age models, considering a Wiener process based degradation. Nevertheless, the estimation procedure we developed highly relies on the non negativity of the gamma process and hence could not be adapted to a general non monotonous Lévy process. The adaptation of the estimation procedure of the present paper to another monotonous Lévy process than the gamma process would be interesting to study.

Hence there remain many estimation procedures to be developed for all these imperfect repair models for deteriorating systems, but the semiparametric estimation of the maintenance eciency for the models mentioned above with various observation schemes, is probably the most challenging problems we aim at investigating in near future.

Extension of the semiparametric estimation method to the Arithmetic Reduction of Degradation model of order innity

The purpose of this section is to prove that the semiparametric estimation method studied previously remains valid in the case where the ARD ∞ model is considered.

First of all, Proposition 6 also holds for the ARD ∞ model, because only the rst two observations are necessary to obtain the model identiability, and both models are the same over [0, 2T ). Now let us recall that the degradation level Y t at any time t such that (j -1)T ≤ t < jT can be expressed as

Y t = X (j) t -X (j) (j-1)T + j-1 k=1 (1 -ρ 0 ) j-k U k where U k = X (k) kT -X (k)
(k-1)T for 1 ≤ k ≤ j -1, and ρ 0 is the true value of the maintenance actions eciency. The observation times are still jT -for 1 ≤ j ≤ n, and let us set Y j = Y jT -. Therefore, we have

Y j = U j + j-1 k=1 (1 -ρ 0 ) j-k U k . (4.20)
Note that the degradation level right after the jth maintenance action is (1 -ρ 0 )Y j , hence

Y j = U j + (1 -ρ 0 )Y j-1
As a result, for each 1 ≤ j ≤ n, the increments U j can be expressed with respect to the observations 4.7. Extension of the semiparametric estimation method to the Arithmetic Reduction of Degradation model of order innity Y j-1 and Y j as well as the maintenance eciency parameter ρ 0 as follows:

U j = Y j -(1 -ρ 0 )Y j-1 .
Similarly as in Subsection 4.3.1, we dene here the function g j as well as the estimator sequence (M j ) 1≤j≤n . It is not necessary to dene the set D j of the possible values for ρ 0 as the sequence (M j ) 1≤j≤n has an explicit expression here.

Let us set g j (ρ, Y ) = Y j -(1 -ρ)Y j-1 , hence we have g j (ρ 0 , Y ) = U j > 0 a.s. Therefore, the parameter ρ 0 fulls the condition g j (ρ 0 , Y ) ≥ 0 for all 1 ≤ j ≤ n. Let us recall that g 1 (ρ 0 , Y ) = Y 1 ≥ 0, and for j = 2, . . . , n the condition

g j (ρ 0 , Y ) ≥ 0 is equivalent to ρ 0 ≥ 1 -Y j /Y j-1 , thus ρ 0 ≥ max 2≤j≤n 1 - Y j Y j-1 .
Finally, we dene the sequence (M j ) 1≤j≤n as

   M 1 = 0 M j = max M j-1 , 1 - Yj Yj-1
for all 2 ≤ j ≤ n which is non-decreasing and bounded by ρ 0 because ρ 0 ∈ [M j , 1) for all 1 ≤ j ≤ n.

The last step is to prove that Lemma 3 still holds here. As a result, all the technical results displayed in Subsection 4.3.2 hold regardless if the order of the ARD model is 1 of ∞, because they only depend on the sequence properties as well as this lemma. Let us now recall the lemma and prove it within the ARD ∞ context. Lemma 3. Let us set ρ 0 ∈ [0, 1), then g j (ρ, Y ) ≥ 0 implies that ρ 0 -ρ ≤ Uj Uj-1 for all 2 ≤ j ≤ n. Proof. Assume j ≥ 2. We derive from Equation (4.20) that

1 - Y j Y j-1 = 1 - U j + (1 -ρ 0 )Y j-1 Y j-1 = 1 - U j Y j-1 -(1 -ρ 0 )
Hence, ρ 0 -ρ ≤ U j /Y j-1 , and replacing Y j-1 by its expression allows to write that

ρ 0 -ρ ≤ U j U j-1 + (1 -ρ 0 )Y j-2 which leads to the result since (1 -ρ 0 )Y j-2 > 0 a.s.
From this, it can be concluded that the main results displayed in Theorems 1, 2, 3 and 4 are veried whether the underlying model is the ARD 1 or the ARD ∞ .

Chapter 5

Introduction

A system is considered whose intrinsic deterioration is modelled by a gamma process (X t ) t≥0 with shape function a and scale parameter b. The system is subject to periodic (period T ) and instantaneous imperfect maintenance actions, and we place ourselves into the framework of a Arithmetic Reduction of Age model of order one (ARA 1 ). Within the context of virtual age models, an imperfect repair puts back the system to a similar state as it was before, namely the system is rejuvenated. In the case of the ARA 1 model, a repair removes a proportion ρ ∈ [0, 1) of the age of the system accumulated since the last maintenance action. This model is dened below.

In Chapter 6, several estimation methods are developed, starting with the most usual ones, namely the methods of moments and the Maximum Likelihood Estimation (MLE). However, the expression of the joint density of the observations, and consequently the likelihood, is a product of integrals of large dimension, and thus numerical estimations becomes dicult to compute in a classical way. At rst, an expectation maximization algorithm, as well as a dierential evolution maximization (see [START_REF] Storn | Dierential Evolution -A Simple and Ecient Heuristic for Global Optimization over Continuous Spaces[END_REF]) associated with quasi Monte Carlo approximations for the integrals, were developed and tested. However, none of these methods were suitable due to numerical issues and high CPU times. Hence, the MLE is treated here by approximating the integrals by the Monte Carlo and randomized Quasi Monte Carlo methods.

Besides these methods, in order to avoid the problem of high-dimensional integrals, alternative maximum likelihood methods are developed: the composite maximum likelihood and the half data method.

The rst method consists in assuming that the observations are independent, while the second one only takes into account one out of two observations. This allows to reduce the dimension of the integrals to 1.

After the methods study, some illustrations of their performances are provided in Chapter 7.

Let us rst dene the ARA 1 model. Let X (j) j∈N * be a sequence of independent copies of the gamma process X = (X t ) t≥0 , which corresponds to the intrinsic degradation of the system. Let (Z t ) t≥0 stands for the eective degradation level of the system. In the framework of the ARA 1 model, each maintenance action remove ρT units of time to the virtual age of the system, which can be expressed as follows:

V (t, ρ) = t -ρnT for nT ≤ t < (n + 1) T
where V (t, ρ) is the virtual age of the system at real time t. Note that V (nT, ρ) = (1 -ρ) nT and V (nT -, ρ) = [(1 -ρ) n + ρ] T . Also, the nth maintenance action puts the system back to the state it was ρT units of time before, that is at real time nT -ρT = (n -ρ) T . The related virtual age is

V ((n -ρ) T, ρ) = (n -ρ) T -ρ (n -1) T = (1 -ρ) nT = V (nT, ρ) .

Chapter 5. Introduction

This is consistent because at real time nT , the system is put back to the exact state it was at real time (n -ρ) T . This leads to the following expression for the degradation level right after the nth maintenance action:

Z nT = Z (n-ρ)T (5.1)
and nally the general expression of the virtual age of the system for all t ≥ 0 is

V (t, ρ) = j≥0 (t -ρjT ) 1 [jT,(j+1)T ) (t).
Hence, the degradation level Z t for nT < t ≤ (n + 1)T , with n in N * , can be expressed as follows

Z t = Z nT + X (n+1) V (t,ρ) -X (n+1) V (nT,ρ) = Z nT + X (n+1) t-ρnT -X (n+1) (1-ρ)nT with Z nT = Z (n-1)T + X (n) V (nT,ρ) -X (n) V ((n-1)T,ρ) = Z (n-1)T + X (n) (1-ρ)nT -X (n) (1-ρ)(n-1)T ,
and where

X (n+1) t-ρnT -X (n+1) (1-ρ)nT is gamma distributed Γ a(t -ρnT ) -a((1 -ρ)nT ), b . Note that Z nT
can be written as

Z nT = n j=1 X (j) (1-ρ)jT -X (j) (1-ρ)(j-1)T ,
which is also gamma distributed with shape and scale parameters a((1 -ρ)nT ) and b respectively. Thus Z t is the sum of two gamma distributed random variables with the same scale parameter, which entails that Z t is gamma distributed Γ a(t -ρnT ), b for all nT ≤ t < (n + 1)T , and its mean and variance are

given by

E(Z t ) = a(t -ρnT ) b and V(Z t ) = a(t -ρnT ) b 2 .
Based on the same framework as in Example 1, Figure 1.3 represents both the degradation evolution mean and variance of the maintained system, as well as those of the intrinsic deterioration mean (unmaintained system). The scale parameter is here equal to one, hence the mean and the variance are equal as well due to their similar expressions.

Chapter 6

Parametric inference for the Arithmetic Reduction of Age model of order one

Preliminary

As mentioned in the introduction of this part, several estimation methods are developed in this chapter, in order to estimate the ARA 1 model parameters. The periodicity T is assumed to be known, thus the parameters of interest are the parameters of the shape function, the scale parameter b and the maintenance eciency ρ. The estimation methods are developed, and then tested in the next chapter, within the framework of a power law shape function. To be more precise, the shape function is dened as a : t → αt β with α, β > 0. We set θ = (α, β, ρ, b) ∈ Θ the parameter set, with the parameter space Θ = (0, ∞) 2 × [0, 1), ×(0, ∞).

We assume that the degradation level is measured right before the rst n maintenance actions for some n in N * , that is at times T -, 2T -,. . . ,nT -. Also, s i.i.d. systems are considered. For sake of readability, let us dene the following notations:

Z j = Z jT -for 1 ≤ j ≤ n; Z = Z jT -1≤j≤n = (Z j ) 1≤j≤n ; z (i) j
is the observed degradation level of the ith maintained system at times jT -for 1 ≤ j ≤ n and 1 ≤ i ≤ s, which is a realisation of the r.v. Z j ;

z (i) = z (i) j 1≤j≤n
is the complete observations set related to the ith maintained system;

z is the complete observations set, that is z = z (i) j 1≤j≤n, 1≤i≤s
; ξ is the parameter set excluding b, that is ξ = (α, β, ρ).

Note that the Z j 's are mutually dependent. Hence, in the next sections, some methods will be developed considering the observations dierence (increments), that is ∆Z = (Z j -Z j-1 ) 1≤j≤n , with Z 0 = 0 by assumption on the initial deterioration level, and these random variables Z j -Z j-1 are only pairwise dependent.

In the following, for sake of simplicity, the same notations are kept (∆ is omitted), that is 6.1. Preliminary

U j = X (j) V ((j-ρ)T,ρ) -X (j)
V ((j-1)T,ρ) .

Thus, we derive that U j is gamma distributed Γ (a [V ((j -ρ) T, ρ)] -a [V ((j -1) T, ρ)] , b), which matches the denition of U j as

a [V ((j -ρ) T, ρ)] -a [V ((j -1) T, ρ)] = a [(1 -ρ) jT ] -a [(1 -ρ) (j -1) T ] = µ(ξ).
Note that from Equation (5.1), we have

U j = Z jT -Z (j-1)T . (6.1)
Now, let V j be the increment of the underlying degradation process over ((j -ρ) T, jT ), with

V j = Z jT --Z (j-ρ)T = X (j) V (jT -,ρ) -X (j) V ((j-ρ)T,ρ)
and so V j is gamma distributed Γ (ν(ξ), b). Once again from Equation (5.1), we have

V j = Z jT --Z jT . (6.2)
Therefore, Equations (6.1) and (6.2), provide the following expression for the observations increments:

Z jT --Z (j-1)T -= Z jT --Z jT + Z jT -Z (j-1)T -Z (j-1)T --Z (j-1)T = U j + V j -V j-1 = W j -V j-1
where W j = U j + V j and where the random variables U 1 , . . . , U n ,V 1 , . . . , V n are independent. Now the model identiability must be veried. To this aim, we rst need the result given in the lemma below.

Lemma 8. Let X 1 , X 2 , Y 1 and Y 2 be independent and gamma distributed random variables with para-

meters (a 1 , b 1 ), (a 2 , b 1 ), (c 1 , b 2 ) and (c 2 , b 2 ) respectively. If X 1 -X 2 D = Y 1 -Y 2 , then          b 1 = b 2 a 1 = c 1 a 2 = c 2 Proof. Assume that X 1 -X 2 D = Y 1 -Y 2 .
We rst look for an expression of the third centered moment of X 1 . Let us recall that the third moment of the dierence of two random variables is equal to the dierence of the moments. Let µ be the Lévy's measure associated to the random variable X 1 . By [11, Eq. 2.32] and [START_REF] Cont | Financial Modelling With Jump Processes[END_REF]Prop. 3.13], we have

E (X t -E (X t )) 3 = R x 3 µ(dx) = ∞ 0 x 3 a 1 e -b1x
x dx

and thus

E (X t -E (X t )) 3 = 2a 1 b 3 1
The expectations, variances and third centered moments of X 1 -X 2 and Y 1 -Y 2 are equal, which leads to

                           a 1 b 1 - a 2 b 1 = c 1 b 2 - c 2 b 2 a 1 b 2 1 + a 2 b 2 1 = c 1 b 2 2 + c 2 b 2 2 2a 1 b 3 1 - 2a 2 b 3 1 = 2c 1 b 3 2 - 2c 2 b 3 2
The rst and third equation provide b 1 = b 2 , and then from the rst and second equation we have

   a 1 = c 1 a 2 = c 2
which nishes the proof.

Theorem 5. Let Z and Z be two random vectors based on the ARA 1 repair model with parameter θ and θ respectively, with θ,

θ ∈ Θ = (0, ∞) 2 × (0, 1) × (0, ∞). If Z D = Z and n ≥ 2, then θ = θ. Proof. Assume that Z D = Z and n = 2. Then W 1 D = W1 and W 2 -V 1 D = W2 -Ṽ1
, and from Lemma 8

we deduce

                 b = b ω 1 (ξ) = ω 1 ( ξ) ν 1 (ξ) = ν 1 ( ξ) ω 2 (ξ) = ω 2 ( ξ)
with ξ = (α, β, ρ) and ξ = (α, β, ρ). The expressions of ω 1 , ω 2 and ν 1 lead to 6.2. Method of moments estimators

                 b = b αT β = αT β 1 -(1 -ρ) β αT β = 1 -(1 -ρ) β αT β (2 -ρ) β -(1 -ρ) β αT β = (2 -ρ) β -(1 -ρ) β αT β
which can be written as

                 αT β = αT β b = b (1 -ρ) β = (1 -ρ) β (2 -ρ) β = (2 -ρ) β
Let us set r = 1 -ρ ∈ (0, 1), the last two equations become: However, the function r → log(1 + r)/ log(r) is injective over (0, 1), hence necessarily r = r. From this we deduce that ρ = ρ, which leads to β = β and then α = α. Therefore θ = θ.

   r β = r β (1 + r) β = (1 + r)
In other words, the model parameters are identiable as soon as the observations are conducted twice, at times T -and 2T -. Now, we go on with the development of parametric estimation methods.

Method of moments estimators

Here we deal with the Moments Estimation (ME) method, which is the exact copy of the method applied to the ARD models. Hence, very few details are given in the following. See Section 2.2 for the detailed denition of the method as well as its application. First, let us recall that this method is based on the minimization of the distance function D dened as Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one

D(θ , θ 0 ) = d i=1 n j=1 m i (θ , jT -) -m i (θ 0 , jT -) 2
and we rst check the parameters identiability in this framework through the following proposition.

Proposition 8. The parameters of the ARA 1 model are identiable from the ME method, that is D(θ , θ 0 ) = 0 implies that θ = θ 0 for all θ, θ 0 ∈ Θ, as soon as d ≥ 2 and n ≥ 3.

Proof. Here we prove that the assertion D(θ , θ 0 ) = 0 ⇒ θ = θ 0 is true ∀ θ, θ 0 ∈ Θ as soon as n ≥ 3 and d ≥ 2. Assume that D(θ , θ 0 ) = 0, and let θ = (α , β , b , ρ) and θ 0 = (α 0 , β 0 , b 0 , ρ 0 ) be in Θ. So given the function D, we have for all j in {1 , ... , n}:

   m 1 (θ , jT -) = m 1 (θ 0 , jT -) m 2 (θ , jT -) = m 2 (θ 0 , jT -)
which is equivalent to

             α (j -ρ(j -1)) β T β b = α 0 (j -ρ 0 (j -1)) β0 T β0 b 0 α (j -ρ(j -1)) β T β b 2 = α 0 (j -ρ 0 (j -1)) β0 T β0 b 2 0 When j = 1, we have        αT β b = α 0 T β0 b 0 αT β b 2 = α 0 T β0 b 2 0 which leads to    αT β = α 0 T β0 b = b 0 (6.3)
Since the expectations and the variances have similar expression, considering both of them only allows the identiability of the scale parameter and the quantity αT β . Hence only the expectation is considered in the following, and taking into account the second and third observation times leads to

             α (2 -ρ) β T β b = α 0 (2 -ρ 0 ) β0 T β0 b 0 α (3 -2ρ) β T β b = α 0 (3 -2ρ 0 ) β0 T β0 b 0
and from System (6.3) we deduce 6.2. Method of moments estimators

   (2 -ρ) β = (2 -ρ 0 ) β0 (3 -2ρ) β = (3 -2ρ 0 ) β0
The parameters ρ, ρ 0 ∈ [0, 1), thus the quantities 3 -2ρ and 3 -2ρ 0 are strictly positive and we have

log(2 -ρ) log(3 -2ρ) = log(2 -ρ 0 ) log(3 -2ρ 0 )
and because the function

f : → log(2 -x) log(3 -2x)
is injective over [0, 1), we conclude that θ = θ 0 .

In other words, the identiability holds from the ME method as soon as the systems are observed three times, at times T -, 2T -and 3T -, and if at least the rst two moments (expectation and variance) are used. Then, an estimation θ of θ 0 can be obtained through

D(θ) = D( θ) = µ 2 + η 2   n j=1 j -(j -1)ρ β T β 2   -2µ   n j=1 j -(j -1)ρ β T β m1 jT -   -η   n j=1 j -(j -1)ρ β T β m2 jT -   + n j=1 m2 1 jT -+ m2 2 jT -
where m1 (jT -) and m2 (jT -) are the empirical expectations and variances of Z jT -. Now we look for zeros of the partial derivatives of D with respect to µ and η:

   ∂ µ D( θ) = 0 ∂ η D( θ) = 0
which leads to Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one

                                       µ = n j=1 j -(j -1)ρ β T β m1 jT - n j=1 j -(j -1)ρ β T β 2 ρ = n j=1 j -(j -1)ρ β T β m2 jT - n j=1 j -(j -1)ρ β T β 2
Given these expressions for µ and η, the function D( θ) can be written as

  n j=1 m2 1 jT -+ m2 2 jT -   -   n j=1 j -(j -1)ρ β T β m1 jT -   2 +   n j=1 j -(j -1)ρ β T β m2 jT -   2 n j=1 j -(j -1)ρ β T β 2
which becomes a function only depending on β and ρ. In order to obtain estimates β and ρ of β and ρ, it remains to minimize numerically this function, which is equivalent to maximize the quantity

  n j=1 j -(j -1)ρ β T β m1 jT -   2 +   n j=1 j -(j -1)ρ β T β m2 jT -   2 n j=1 j -(j -1)ρ β T β 2 (6.4)
Finally, the estimator θ of θ has the following expression

           n j=1 j -(j -1)ρ β T β m1 jT -   2   n j=1 j -(j -1)ρ β T β m2 jT -     n j=1 j -(j -1)ρ β T β 2   , β , ρ , n j=1 j -(j -1)ρ β T β m1 jT - n j=1 j -(j -1)ρ β T β m2 jT -         

Maximum Likelihood Estimation

In this section we deal with the MLE method. We begin with the case where one single system is observed (s = 1) and we start by computing the joint p.d.f. of the random vector Z, from where the likelihood function is easily derived when s > 1.

Maximum Likelihood Estimation

We know from Theorem 5 that the model identiability holds only if n ≥ 2. Then, let us set n ≥ 2, the joint p.d.f. f Z of Z is given by

f Z (Z) = R n-1 f Z|V n-1 (Z | v 1 , . . . , v n-1 ) f V n-1 (v 1 , . . . , v n-1 ) dv 1 . . . dv n-1 where Z = (z 1 , z 2 , • • • , z n ) ∈ R + × R n-1 . The r.v. V j are independent, then f V n-1 (v 1 , . . . , v n-1 ) = n-1 j=1 f Vj (v j )
and because Z j = U j + V j -V j-1 , we can express the p.d.f. f Z|V n-1 as

f Z1|V n-1 (z 1 | v 1 , . . . , v n-1 ) = f U1+V1|V1 (z 1 | v 1 ) = f U1 (z 1 -v 1 ) , f Zj |V n-1 (z j | v 1 , . . . , v n-1 ) = f Uj +Vj -Vj-1|Vj-1,Vj (z j | v j-1 , v j ) = f Uj (z j + v j-1 -v j )
for 1 < j < n, and

f Zn|V n-1 (z n | v 1 , . . . , v n-1 ) = f Un+Vn-Vn-1|Vn-1 (z n | v n-1 ) = f Wn (z n + v n-1 ) .
Therefore, with the convention v 0 = 0, the joint p.d.f. f Z of Z can be expressed as follows:

f Z (Z) = R n-1   n-1 j=1 f Vj (v j ) f Uj (z j + v j-1 -v j )   × f Wn (z n + v n-1 ) dv 1 . . . dv n-1 . (6.5)
We now reduce the integration domain to [0, 1] n-1 in order to be able to apply either Monte Carlo (MC) or Quasi Monte Carlo (QMC) methods for the integrals computations. This will prevent issues within the implementation of this method, since the computation of an integral over R n-1 might be an issue when n is large.

Proposition 9. The joint p.d.f. of the random vector Z can be expressed as

f Z (Z) = n-1 j=1 C j -m j+1 [0,1] n-1 G θ (x, z) dx 1 . . . dx n-1
where x = (x 1 . . . x n-1 ), m j = max(0, -z j ), C j = j i=1 z i , and

G θ (x, z) = n-1 j=1 f Vj ((C j -m j+1 ) x j + m j+1 )
Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one

× n-1 j=1 f Uj (z j + (C j-1 -m j ) x j-1 -(C j -m j+1 ) x j + m j -m j+1 ) × f Wn z n + x n-1 C n-1 -m n + m n
Proof. The proof stands in two steps. First, we prove that the integral domain can be reduced to the set n-1 j=1 max(0, -z j+1 ) , C j for all 1 ≤ j ≤ n -1. Then, a variable change will ensure the result.

The function to integrate in Equation (6.5) is a product of p.d.f., which can be zero over R n-1 . Hence, in order to reduce the integral domain, we look for a subset of R n-1 containing all the value of the v j for which the p.d.f. are non zero. Here, the functions f Uj , f Vj and f Wn are all non zero if

0 < v j for 1 ≤ j ≤ n -1; 0 < z j + v j-1 -v j for 1 ≤ j ≤ n -1; max(0, -z n ) < v n-1 .
Let us set m j = max(0, -z j ) and C j = j i=1 z i > 0 for 1 ≤ j ≤ n -1, and show by induction that

m j+1 < v j < C j for all 1 ≤ j ≤ n -1. If j = 1, then 0 < z 1 + v 0 -v 1 and v 1 < z 1 = C 1 . Similarly, 0 < z 2 + v 1 -v 2 hence v 2 -z 2 < v 1 , and because v 2 > 0 we have -z 2 < v 1 , and 0 < v 1 implies that m 2 < v 1 . Hence m 2 < v 1 < C 1 .
Now assume that m j < v j-1 < C j-1 for 2 ≤ j ≤ n -2, the aim is to show that m j+1 < v j < C j . We know by assumption that

v j , v j+1 > 0 0 < z j + v j-1 -v j 0 < z j+1 + v j -v j+1
and v j < z j + C j-1 = C j because v j < z j + v j-1 . Also, v j+1 -z j+1 < v j , v j > 0 and v j-1 > 0, which implies that -z j+1 < v j and nally m j+1 < v j because of the positivity of v j .

Finally, if j = n -1 we know that m n < v n-1 . Moreover, v n-1 < z n-1 + v n-2 < z n-1 + C n-2 = C n-1 , which leads to m n < v n-1 < C n-1 .
We conclude that the integral domain can be written as

n-1 j=1 max(0, -z j+1 ) , C j .
Note that the p.d.f. might be zero within this domain, but cannot be non zero outside. Now, the following change of variables can be applied. Let us set

x j = v j -m j+1 C j -m j+1 for all 1 ≤ j ≤ n -1 hence dv j = (C j -m j+1 ) dx j
6.4. Half data method for 1 ≤ j ≤ n -1, and the integration domain becomes [0, 1] n-1 . Finally, substituting the v j and the dv j by x j and dx j leads to the result. We now consider s observed systems. Let us set S (z i,n ) such that

S (z i,n ) = n-1 j=1 C i,j -m i,j+1
with C i,j = j k=1 z i,k and m i,j+1 = max (0, -z i,j+1 ). The log-likelihood is easily derived from the previous result, which has the following expression

(θ | z) = s i=1 log S (z i,n ) + s i=1 log [0,1] n-1 G θ (x, z) dx 1 . . . dx n-1 .
(6.6)

Half data method

We recall that the random variables Z j = U j + V j -V j-1 for 1 ≤ j ≤ n are not independent, which leads to a cumbersome expression for the joint p.d.f., which is dicult to assess from a numerical point of view, as we have seen in the previous section. However, the Z j 's with non consecutive indices j's can be seen to be independent, as these variables are pairwise dependent. We hence follow the idea of the sub-sample methods (see [START_REF] Babu | Subsample and half-sample methods[END_REF][START_REF] Hartigan | Using subsample values as typical values[END_REF]), and we consider only one-out-of two Z j 's, restricting the information to either odd indices, or even indices, and letting down the other ones.

Half data based on the odd indexes

For sake of simplicity, n is assumed here to be odd, and we set p = (n + 1)/2. The half data vector Z is composed by the r.v.s. Z j whose index is odd, that is

Z = Z 2j-1 1≤j≤p = W 2j-1 -V 2j-2 1≤j≤p
where the Z 2j-1 are independent. Based on Lemma 8, the identiability holds if the following system of equations

         ω 2j-1 (ξ) = ω 2j-1 (ξ 0 ) ν 2j-2 (ξ) = ν 2j-2 (ξ 0 ) (6.7)
for j = 1, . . . , p has a unique solution in (R + ) 2 × (0, 1), that is ξ = ξ 0 . We rst consider the case n = 3

(where two observations per trajectory are used), which will be seen to be non identiable on the whole domain. We next go the the case n=5 (where three observations per trajectory are used), which will be shown to provide identiability on the whole domain.

Case n = 3 System (6.7) becomes :

                       αT β = α 0 T β0 (3 -2ρ) β -(2 -2ρ) β αT β = (3 -2ρ 0 ) β0 -(2 -2ρ 0 ) β0 α 0 T β0 (2 -ρ) β -(2 -2ρ) β αT β = (2 -ρ 0 ) β0 -(2 -2ρ 0 ) β0 α 0 T β0
which is equivalent to

                       αT β = α 0 T β0 (1 + 2r) β -(2r) β = (1 + 2r 0 ) β0 -(2r 0 ) β0 (1 + r) β -(2r) β = (1 + r 0 ) β0 -(2r 0 ) β0
with r = 1 -ρ and r 0 = 1 -ρ 0 . We concentrate on the last two equations and try to identify the cases where they imply that (r, β) = (r 0 , β 0 ). In such cases, the rst equation leads to α = α 0 , and identiability is obtained.

Let us set

h 1 (r, β) = (1 + 2r) β -(2r) β h 2 (r, β) = (1 + r) β -(2r)
β with β > 0 and r ∈ (0, 1).

Lemma 9. Let (r, β) and (r 0 , β 0 ) be such that h i (r, β) = h i (r 0 , β 0 ) for i = 1, 2. Then,

1. if β 0 = β or r = r 0 then (r, β) = (r 0 , β 0 ); 2. if β 0 ≤ 1 then (r, β) = (r 0 , β 0 ); 3. if β 0 > 1 then β > 1; 4. if β 0 , β > 1, then β 0 > β implies r > r 0 and β 0 < β implies r < r 0 ; 5. if (r 0 , β 0 ) is such that h 2 (r 0 , β 0 ) < 1, then (r, β) = (r 0 , β 0 ).
Proof. The partial derivatives of h 1 and h 2 are given by

∂ r h 1 (r, β) = 2β (2r + 1) β-1 -(2r) β-1 ∂ β h 1 (r, β) = (2r + 1) β ln (2r + 1) -(2r) β ln (2r) ∂ r h 2 (r, β) = β (r + 1) β-1 -2 (2r) β-1 ∂ β h 2 (r, β) = (r + 1)
β ln (r + 1) -(2r) β ln (2r)

Half data method

Note that both partial derivatives with respect to β are strictly positive. Moreover, ∂ r h 1 (r, 1) = 0 and the quantity ∂ r h 1 (r, β) has the same sign as β -1, and if β ≤ 1 then ∂ r h 2 (r, β) < 0.

Also, when β > 1, the partial derivative ∂ r h 2 (r, β) changes sign at g(β) from positive to negative, where

g(β) = 1 2 β β-1 -1
and β → g(β) is non-decreasing from (1, ∞) to (0, 1).

1. Assume β = β 0 . If β ≤ 1, then r -→ h 2 (r, β) is injective and thus h 2 (r, β) = h 2 (r 0 , β 0 ) implies r = r 0 . Similarly, if β > 1 then r -→ h 1 (r, β) is injective and r = r 0 . Assume now that r = r 0 , hence β = β 0 because the function β -→ h 1 (r, β) is injective.
2. Assume β 0 = 1, thus h 1 (r 0 , 1) = 1 and by assumption h 1 (r, β) = 1. However if β > 1, then h 1 (r, β) > h 1 (0, β) = 1 and in the same way, if β < 1, then h 1 (r, β) < 1. This implies β = β 0 = 1, and r = r 0 is directly obtained from the previous case.

3. Assume β 0 < 1, thus h 1 (r 0 , β 0 ) < h 1 (r 0 , 1) = 1 = h 1 (r, 1) and so h 1 (r, β) < 1, hence β < 1.

If β > β 0 , then h 1 (r, β) > h 1 (r, β 0 ). Moreover h 1 (r, β) = h 1 (r 0 , β 0 ) hence 0 = h 1 (r, β) -h 1 (r 0 , β 0 ) = h 1 (r, β) -h 1 (r, β 0 ) + h 1 (r, β 0 ) -h 1 (r 0 , β 0 )
and so h 1 (r, β 0 ) -h 1 (r 0 , β 0 ) < 0 because h 1 (r, β) > h 1 (r, β 0 ), which leads to r > r 0 .

Let us now dene h 3 such that h 3 (r, β) = h 1 (r, β)-h 2 (r, β). We have h 3 (r, β) = (1+2r) β -(1+r) β and

∂ r h 3 (r, β) = β 2 (2r + 1) β-1 -(1 + r) β-1 > 0; ∂ β h 3 (r, β) = (2r + 1) β ln (2r + 1) -(1 + r) β ln (1 + r) > 0.
Based on the fact that β > β 0 and r > r 0 , it can be stated that

h 3 (r, β) > h 3 (r 0 , β 0 ) ,
which is impossible because h 3 (r, β) = h 3 (r 0 , β 0 ) by denition, hence (r, β) = (r 0 , β 0 ). In the same way, this result holds when β < β 0 .

Assume β

0 > 1, then h 1 (r, β) = h 1 (r 0 , β 0 ) > h 1 (r 0 , 1) = 1 = h 1 (r, 1) hence β > 1.
With similar arguments as in the previous case, if β < β 0 , then r > r 0 . Finally if β > β 0 , then r < r 0 by symmetry.

5. If β 0 ≤ 1 then (r, β) = (r 0 , β 0 ) from the second case.

Assume β 0 > 1 and h 2 (r 0 , β 0 ) < 1. Then, we know from the third case that β > 1, which implies that r → h 2 (r, β) is increasing over (0, g(β)]. Also, note that h 2 (0, β) = 1 for all β > 0, and Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one because r → h 2 (r, β) is increasing over (0, g(β)], necessarily h 2 (r, β) > 1, which is impossible as h 2 (r 0 , β 0 ) = h 2 (r, β). Hence r 0 > g(β 0 ) and r > g(β) with similar arguments. Also, either β 0 > β and r > r 0 , or β 0 < β and r < r 0 . If β 0 > β and r > r 0 then h 2 (r 0 , β 0 ) > h 2 (r 0 , β) because ∂ β h 2 (r, β) > 0. Moreover, r > r 0 > g(β 0 ) > g(β) and the function r → h 2 (r, β) is decreasing over [g(β), 1), hence h 2 (r 0 , β) > h 2 (r, β). This leads to h 2 (r 0 , β 0 ) > h 2 (r, β) which is impossible since h 2 (r 0 , β 0 ) = h 2 (r, β). In the same way, the result holds when β 0 < β and r < r 0 .

Corollary 2. Let (α, β, r, α 0 , β 0 , r 0 ) be such that (6.7) is satised. If β 0 ≤ 1 or h 2 (r 0 , β 0 ) < 1 then ξ = ξ 0 as soon as n ≥ 3.

Proof. From Lemma 9, if β 0 ≤ 1 or h 2 (r 0 , β 0 ) < 1 then (β, r) = (β 0 , r 0 ), and the rst line of (6.7)

ensures that α = α 0 . Let us recall that r = 1 -ρ and r 0 = 1 -ρ 0 , this leads to the conclusion ξ = ξ 0 .

The remaining cases now are 1 < β < β 0 and r < r 0 , or 1 < β 0 < β and r > r 0 .

Let (r, β) and (r 0 , β 0 ) be such that:

(*)

• h i (r, β) = h i (r 0 , β 0 ) for i = 1, 2,
• 1 < β < β 0 and r 0 < r.

Note that the case where 1 < β 0 < β and r > r 0 is identical to this one by symmetry, hence the study of one of the cases is enough. Now, we investigate whether Conditions (*) are compatible. To this end, let us now dene the function

u 2 by u 2 (γ) = (1 + γr) β -(2r) β -(1 + γr 0 ) β0 + (2r 0 ) β0
for γ > 0, where the parameters (r, β, r 0 , β 0 ) are xed and omitted in the notation u 2 (γ), for sake of simplication. From Conditions (*) we have u 2 (1) = u 2 (2) = 0, which is possible only if u 2 is nonmonotonic over [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF][START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. Therefore, in the following, the monotonicity of u 2 is studied.

Lemma 10. Suppose Conditions (*) to hold. Then u 2 (γ) > 0 if and only if rβ r0β0 > w (γ), with

w (γ) = (1 + γr 0 ) β0-1 (1 + γr) β-1
for all γ > 0.

Moreover, w (γ) reaches its minimum at γ 0 given by

γ 0 = r 0 -r -r 0 β 0 + rβ rr 0 (β 0 -β) .
Proof. The derivative of u 2 is given by

u 2 (γ) = rβ (1 + γr) β-1 -r 0 β 0 (1 + γr 0 ) β0-1 ,
thus u 2 (γ) > 0 if and only if rβ r0β0 > w (γ). Then, w (γ) has the same sign as the quantity

(β 0 -1) r 0 (1 + γr) -(β -1) r (1 + γr 0 )
6.4. Half data method which is positive if and only if γ > γ 0 , and so w (γ) > w (γ 0 ) for all γ = γ 0 .

Corollary 3. Suppose Conditions (*) to hold. If the function u 2 is non-monotonic over [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF][START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF], then one of the following assertions holds : Proof. Remember from Lemma 10 that the function w is rst decreasing on [0, γ 0 ] and next increasing on [γ 0 , ∞), and note that there are three possible cases, according to whether

γ 0 ≤ 1, γ 0 ≥ 2 or γ 0 ∈ (1, 2).
Assume rst that γ 0 ≤ 1. Then w is increasing on [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF][START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF], and hence w (1) < w (2). Based on Lemma 10 again, if the function u 2 is non monotonic, it means that rβ r0β0 must belong to (w (1) , w (2)), otherwise, the sign of u 2 would remain a constant on [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF][START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF].

In the same way, if γ 0 ≥ 2, then w is decreasing on [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF][START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF], and hence w (2) < w (1). Also, rβ r0β0 must belong to (w (2) , w (1)), with similar arguments.

Finally, assume that γ 0 ∈ (1, 2). Then, the function w is decreasing on [1, γ 0 ] and increasing on [γ 0 , 2],

so that w ([1, 2]) = [w (γ 0 ) , max (w (1) , w (2))] .
If the function u 2 is non monotonic, it now means that rβ r0β0 must belong to (w (γ 0 ) , max (w (1) , w (2))).

Lemma 11. Suppose Conditions (*) to hold. Then γ 0 > 1 implies rβ r0β0 ≥ w (1).

Proof. Assume that γ 0 > 1. Our aim is to show that rβ r0β0 ≥ w (1) which is equivalent to prove that

g (r, β) ≤ 0 (6.8) with g (r, β) = r 0 β 0 (1 + r 0 ) β0-1 -rβ (1 + r) β-1
, for all r in (0, 1) and β > 0 such that γ 0 > 1.

As a rst step, let us show that γ 0 > 1 implies that

r > v 1 (β) = r 0 β 0 -1 β -1 -r 0 (β 0 -β) •
Let us note that γ 0 > 1 can be written as

r 1 + r > r 0 1 + r 0 β 0 -1 β -1 , (6.9)
Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one or equivalently as

r 1 - r 0 1 + r 0 β 0 -1 β -1 > 1. (6.10)
Let us check that the factor between brackets is positive in the previous expression. As r/ (r + 1) < 1, Inequality (6.9) entails that

1 > r 0 1 + r 0 β 0 -1 β -1 (6.11)
which can be written as

β > β 1 ≡ 1 + r 0 1 + r 0 (β 0 -1) .
Now it is easy to check that based on β > β 1 that Inequality (6.11) is true, and we now obtain from (6.10) that

r > 1 1 -r0 1+r0 β0-1 β-1 -1 = v 1 (β) .
Based on the fact that g is decreasing with respect to r, we now obtain that

g (r, β) < g (v 1 (β) , β) .
In order to show (6.8), it is now sucient to prove that g (v 1 (β) , β) < 0. We have

g (v 1 (β) , β) = r 0 β 0 (1 + r 0 ) β0-1 - r 0 (β 0 -1) β -1 -r 0 (β 0 -β) β 1 + r 0 (β 0 -1) β -1 -r 0 (β 0 -β) β-1 = r 0 β 0 (1 + r 0 ) β0-1 - r 0 (β 0 -1) β -1 -r 0 (β 0 -β) β (r 0 + 1) β-1 (β -1) β-1 (β + βr 0 -r 0 β 0 -1) β-1 = r 0 (r 0 + 1) β-1 β 0 (1 + r 0 ) β0-β - β 0 -1 β -1 β (β -1) β (β + βr 0 -r 0 β 0 -1) β = r 0 (r 0 + 1) β-1 φ (β) with φ (β) = β 0 (1 + r 0 ) β0-β -(β 0 -1) β β -1 1 1 -r 0 β0-β β-1 β .
The quantity g (v 1 (β) , β) has the same sign as φ (β), hence it remains to prove that φ (β) ≤ 0. First, we have β/(β -1) > β 0 /(β 0 -1) because β < β 0 . This implies

φ (β) < β 0   (1 + r 0 ) β0-β - 1 1 -r 0 β0-β β-1 β   
6.4. Half data method and consequently it is sucient to show that

β ln 1 -r 0 β0-β β-1 1 + r 0 < -β 0 ln (1 + r 0 ) . (6.12)
Let us set z such that

z = β 0 -β β -1 •
Inequation (6.12) is equivalent to τ (z) < -β 0 ln (1 + r 0 ) where

τ (z) = z + β 0 z + 1 ln 1 -r 0 z 1 + r 0 with 0 < z < β 0 -β 1 β 1 -1 = 1 r 0 .
The derivative of τ is given by

τ (z) = - β 0 -1 (z + 1) 2 ln 1 -r 0 z 1 + r 0 + z + β 0 z + 1 -r 0 1 -r 0 z
and has the same sign as

σ (z) = -(β 0 -1) ln 1 -r 0 z 1 + r 0 -r 0 (z + β 0 ) (z + 1) 1 -r 0 z . Also σ (z) = -(β 0 -1) -r 0 1 -r 0 z -r 0 1 (1 -r 0 z) 2 -r 0 z 2 + 2z + β 0 + r 0 β 0 + 1 = r 0 (z + 1) zr 0 -r 0 β 0 -2 (zr 0 -1)
2 and σ (z) < 0 because zr 0 < 1. Thus the function σ is non-increasing, with σ (0) < 0, therefore σ (z) < 0. This leads to τ (z) < 0 for all 0 < z < 1 r0 , which implies that τ is also non-increasing, with τ (0) = -β 0 ln (1 + r 0 ). Therefore, Inequation (6.12) is satised because τ (z) < -β 0 ln (1 + r 0 ) for all 0 < z < 1 r0 , which implies g(r, β) ≤ 0 and nally rβ r0β0 ≥ w (1). This is inconsistent with the results of Lemma 11 because γ 0 > 1 implies rβ r0β0 ≥ w(1), hence this point cannot be veried. Moreover, in the case where γ 0 ∈ (1, 2), we have rβ r0β0 ∈ (w (γ 0 ) , max (w(1), w (2))), that is either rβ r0β0 ∈ (w (γ 0 ) , w (2)) or rβ r0β0 ∈ (w (γ 0 ) , w (1)). Once again because γ 0 > 1 implies rβ r0β0 ≥ w(1), the second case is impossible. It remains the case rβ r0β0 ∈ (w (γ 0 ) , w (2)), but w(1) > w (γ 0 ) by denition of γ 0 , hence necessarily rβ r0β0 ∈ (w (1) , w (2)) with w (1) < w (2).

In summary, if h i (r, β) = h i (r 0 , β 0 ) for i = 1, 2, then the identiability holds when β 0 ≤ 1 or h 2 (r 0 , β 0 ) < 1, and in the case where β 0 > 1 we know that β > 1. Moreover, if β < β 0 and r > r 0 we have from Corollary 4: Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one γ 0 < 2; rβ r0β0 ∈ (w (1) , w (2)), with w (1) < w (2).

We now look at the possible zone for (r, β) with a xed (r 0 , β 0 ) based on the previous results, from a numerical point of view.

Example 7. In Figure 6.1 are illustrated the possible zones for (r, β) based on Corollary 4 with (r 0 , β 0 ) = (0.01, 5) and (r 0 , β 0 ) = (0.05, 4). The blue hatched area corresponds to the (r, β) such that rβ r0β0 > w(1), while the grey hatched area corresponds to the (r, β) such that rβ r0β0 < w(2), and in both zones γ 0 < 2. Therefore, from Corollary 4, if there is a solution (r, β) satisfying Conditions (*) then (r, β) belongs to the intersection of the hatched zones, which can be very small as seen in the right-hand plot. > w(1), while the grey hatched area corresponds to the (r, β) such that γ 0 < 2 and rβ r0β0 < w(2). The red crosses represent the values of (r 0 , β 0 ), with (r 0 , β 0 ) = (0.01, 5) on the left-hand plot and (r 0 , β 0 ) = (0.05, 4) on the right-hand one.

Up to now, we only know that, if h i (r, β) = h i (r 0 , β 0 ) for i = 1, 2 with 1 < β < β 0 and r 0 < r, then (r, β) should be in the intersection of the two hatched zones in Figure 6.1, obtained through Corollary 4. Based on some numerical investigations, it seems that such an intersection is never empty whatever (r 0 , β 0 ) is. Hence, we are not able to conclude about possible identiability, based on the previous theoretical results.

We next investigate identiability from a numerical point of view. For each (r 0 , β 0 ) in (0, 1) × (1, ∞), let us dene the function

g (r0,β0) (r, β) = h 1 (r, β) -h 1 (r 0 , β 0 ) 2 + h 2 (r, β) -h 2 (r 0 , β 0 ) 2 ,

Half data method

which is zero as soon as h i (r, β) = h i (r 0 , β 0 ) for i = 1, 2. We begin with a rst numerical illustration of the behavior of this function.

Example 8. Let us set (r 0 , β 0 ) ≈ (0.07, 4.14). The function g (r0,β0) (r, β) is plotted in Figure 6.2. It can be seen that besides (r 0 , β 0 ) where the function g (r0,β0) (r, β) is zero, there is another point (r, β) = (0.1, 3) for which this function is very close to zero, and there is one single such point. From a numerical point of view, there hence seems to exist a solution to Conditions (*), and identiability seems not to hold. The possible area for (r, β) given by Corollary 4 is plotted in Figure 6.3, where (r, β) = (0.1, 3) is indicated by a red cross. We observe that (r, β) = (0.1, 3) is in the possible area that was theoretically obtained, which is coherent. More specically, we have:

h i (r, β) ≈ h i (r 0 , β 0 ) for i = 1, 2;
1 < β < β 0 and r 0 < r; γ 0 ≈ -2.48 < 2; Similar results were found dealing with both points from Example 7. In each case, a unique (r, β) was found: (r, β) ≈ (0.0213, 2.385) when (r 0 , β 0 ) = (0.01, 5) and (r, β) ≈ (0.072, 2.853) when (r 0 , β 0 ) = (0.05, 4).

These numerical counter-examples lead us to conjecture that the identiability does not hold over the entire set (0, 1) × (0, ∞) when n = 3. For some (r 0 , β 0 ) ∈ (0, 1) × (1, 50), approximate solutions (r, β) such that g (r0,β0) (r, β) = 0 are searched.

To be more specic, for each (r 0 , β 0 ), we look for (r, β) such that

), with w (1) < w (2). The red crosses represent (r, β) and (r 0 , β 0 ).

(r, β) = (r 0 , β 0 ); (r, β) is not in a neighbourhood of (r 0 , β 0 ), that is (r 0 , β 0 ) -(r, β) 1 > 0.01; g (r0,β0) (r, β) < 10 -9 .

For each (r 0 , β 0 ) that has been considered in the study, either there was no solution for (r, β), so that it seems that such a (r 0 , β 0 ) belongs to the identiability zone, either there was one single solution (r, β) (as in the previous examples), and identiability does not seem to hold for such a (r 0 , β 0 ). The points where identiability does not seem to hold are indicated in Figure 6.4 with a black cross. The red line corresponds to the boundary of the zone containing such non identiable points, and no solution has been found outside this area. The area where the identiability were proved to hold is also plotted in blue. Besides this zone, the previous numerical results lead us to the conjecture that the identiability does not hold within the red area, while it holds outside. The black points are the values for (r 0 , β 0 ) for which there exists (r, β) such that g (r0,β0) (r, β) < 10 -9 and (r, β) is not in a neighbourhood of (r 0 , β 0 ). The red lines correspond to the boundary of such points, and the area where theoretical identiability holds is plotted in blue.

Case n = 5

The system (6.7) becomes :

                                                 αT β = α 0 T β0 h 1 (r, β) = h 1 (r 0 , β 0 ) h 2 (r, β) = h 2 (r 0 , β 0 ) h 3 (r, β) = h 3 (r 0 , β 0 ) h 4 (r, β) = h 4 (r 0 , β 0 ) (6.13)
where

h 3 (r, β) = (1 + 4r) β -(4r) β h 4 (r, β) = (1 + 3r) β -(4r) β
Here Conditions (*) have to be redened as two equations were added to the previous system.

Let (r, β) and (r 0 , β 0 ) be such that:

(**)

• h i (r, β) = h i (r 0 , β 0 ) for i ∈ {1, 2, 3, 4}
• 1 < β < β 0 and r 0 < r.

As before for n = 3, let us dene u 4 such that

u 4 (γ) = (1 + γr) β -(4r) β -(1 + γr 0 ) β0 + (4r 0 ) β0
for γ ∈ [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF][START_REF] Babu | Subsample and half-sample methods[END_REF]. Hence from Conditions (**) we have

u 2 (1) = u 2 (2) = 0 u 4 (3) = u 4 (4) = 0
Using similar arguments as those in Corollary 3 leads to the following corollary for u 4 .

Corollary 5. Suppose Conditions (**) to hold. If the function u 4 is non-monotonic over [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF][START_REF] Babu | Subsample and half-sample methods[END_REF], then one of the following assertions holds :

γ 0 ≤ 3 and rβ r0β0 ∈ (w (3) , w (4))
;

γ 0 ≥ 4 and rβ r0β0 ∈ (w (4) , w (3));
γ 0 ∈ (3, 4) and rβ r0β0 ∈ (w (γ 0 ) , max (w (3) , w (4))).

In Figure 6.5, the same area as in Figure 6.3 is plotted in blue, as well as the area where the results of Corollary 5 are satised, for (r 0 , β 0 ) ≈ (0.07, 4.14). The intersection of the two areas seems empty, which indicates that identiability holds in this case. The red cross represents (r 0 , β 0 ). Theorem 6. (Identiability) Let (ξ, ξ 0 ) be such that (6.13) is satised, then ξ = ξ 0 as soon as n ≥ 5.

Half data method

Proof. From Lemma 2, if β 0 ≤ 1 then the identiability holds. Now assume Conditions (**) to hold.

Necessarily γ 0 < 2 and rβ r0β0 ∈ (w (1) , w (2)) based on Corollary 4.

If γ 0 ≤ 1, then rβ r0β0 ∈ (w (1) , w (2)) and rβ r0β0 ∈ (w (3) , w (4)) from Corollary 5. However, w (γ) is non-decreasing over [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF][START_REF] Babu | Subsample and half-sample methods[END_REF], which is inconsistent. Hence γ 0 ∈ (1, 2).

Then, rβ r0β0 ∈ (w (1) , w (2)) and rβ r0β0 ∈ (w (3) , w (4)) and thus rβ r0β0 ∈ (max (w (1) , w (3)) , min (w (2) , w (4))). The function w is non-decreasing over (γ 0 , 4), which leads to w (4) > w (3) > w (2) as γ 0 < 2, and we know that w (1) < w (2). Hence max (w (1) , w (3)) = w (3) and min (w (2) , w (4)) = w (2). Therefore rβ r0β0 ∈ (w (3) , w (2)), which is impossible because w (3) > w (2). As a conclusion, the case where 1 < β < β 0 and r > r 0 is impossible, as well as the case β > β 0 > 1 and r < r 0 by symmetry, thus (β, r) = (β 0 , r 0 ) and nally ξ = ξ 0 .

Likelihood expression

By denition, these random variables Z 2j-1 , for 1 ≤ j ≤ p, are independent. Hence the joint density function of Z is given by

f Z (z n ) = p j=1 f Z2j-1 (z 2j-1 ) = f W1 (z 1 ) p j=2 f W2j-1-V2j-2 (z 2j-1 ). (6.14)
For 2 ≤ j ≤ p, the p.d.f. of W 2j-1 -V 2j-2 can be developed as

f W2j-1-V2j-2 (z 2j-1 ) = R f W2j-1 (z 2j-1 + x)f V2j-2 (x) dx.
Finally, replacing the p.d.f.'s by their expressions in Equation (6.14), we obtain the following expression for the joint density of Z

f Z (z) = z ω1-1 1 b Ωp exp   -b p j=1 z 2j-1   Γ(ω 1 ) p j=2 Γ(ω 2j-1 )Γ(ν 2j-2 ) ×   p j=2 ∞ max(0,-z2j-1) x ν2j-2-1 (z 2j-1 + x) ω2j-1-1 exp (-2bx) dx   1 R + (z 1 )
with ν 0 = 0 and Ω p = p j=1 ω 2j-1 + ν 2j-2 . Now s identical systems are considered. Hence an index i is added to each quantity, referring to the ith system, and the observation set is given by

Z = Z (i) 2j-1 1≤i≤s , 1≤j≤p
Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one whose related density function is given by

f Z (z) = b s Ωp exp   -b s i=1 p j=1 z (i) 2j-1   s i=1 z (i) 1 ω1-1 Γ(ω 1 ) s p j=2 Γ(ω 2j-1 ) s Γ(ν 2j-2 ) s ×   s i=1 p j=2 ∞ max(0,-z (i) 2j-1 )
x ν2j-2-1 (z

(i) 2j-1 + x) ω2j-1-1 exp (-2bx) dx   1 R + (z 1 ).
Finally, the log-likelihood has the following expression

(θ | z) = s Ωp log b -b s i=1 p j=1 z (i) 2j-1 + (ω 1 -1) s i=1 log z (i) 1 -s log Γ(ω 1 ) (6.15)
-

s p j=2 log Γ(ω 2j-1 ) -s p j=2 log Γ(ν 2j-2 ) + s i=1 p j=2 log I θ | z (i) 2j-1
where

I θ | z (i) 2j-1 = ∞ max (0,-z (i) 2j-1 )
x ν2j-2-1 (z

(i) 2j-1 + x) ω2j-1-1 exp (-2bx) dx.
Note that all the results hold if n is an even integer, but in this case p must be such that p = n/2.

Half data based on the even indexes

In a manner similar to that used in the previous section, we only consider the Z j 's whose index is even, in order for these random variables to be independent. The identiability is not studied here. Indeed, in Chapter 7 dealing with numerical analysis, we will see that this method may have identiability issues.

Here n is assumed to be even, and we set p = n/2. The considered observations set is

Z = (Z j ) 1≤j≤p = (W j -V j-1 ) 1≤j≤p
The joint p.d.f. of Z is here given by

f Z (z n ) = p j=1 f Z2j (z 2j ) = p j=1 f W2j -V2j-1 (z 2j )
Therefore, in a similar way than for the previous method, we derive the joint density expression below

f Z (z) = b Ω p exp   -b p j=1 z 2j   p j=1 Γ(ω 2j )Γ(ν 2j-1 )   p j=1 ∞ max(0,-z2j ) x ν2j-1-1 (z 2j + x) ω2j -1 exp (-2bx) dx   (6.16)
where 6.5. Maximum composite likelihood estimation

Ω p = p j=1 ω 2j + ν 2j-1 .
Now, considering s i.i.d. systems, we easily derive the log-likelihood expression from Equations (6.15) and (6.16), which is given by

(θ | z) = s Ωp log b -b s i=1 p j=1 z (i) 2j -s p j=1 log Γ(ω 2j ) (6.17) -s p j=1 log Γ(ν 2j-1 ) + s i=1 p j=1 log I θ | z (i) 2j
where

I (θ | z i,2j ) = ∞ max(0,-z (i) 2j )
x ν2j-1-1 (z

(i) 2j-1 + x) ω2j -1 exp (-2bx) dx.
Note that if n is odd, the previous results hold with p = (n -1)/2.

Maximum composite likelihood estimation

Another way to avoid dependency issues is to assume that the variables are independent, even if they are not. It is the purpose in this section, where two MLE methods are developed considering composite likelihood. To be more precise, the rst one deal with the observations (assumed mutually independent), while the second one is based on the increments (also assumed mutually independent), according to the denition provided by [START_REF] Varin | Composite likelihood methods[END_REF].

Composite likelihood based on the observations

For all j ∈ {1, . . . , n}, the r.v. Z j = Z jT -are gamma distributed Γ a (jT -ρ(j -1)T ) , b . Let us set a j = a (jT -ρ(j -1)T ) the shape parameters of these distributions, we have the following expression for the p.d.f. of

Z j f Zj (z) = b aj Γ(a j ) z aj -1 exp (-bz) 1 R + (z)
Assuming that these r.v.s are independent, the likelihood is given by

L (z | θ) = b n j=1 aj n j=1 Γ(a j ) exp   -b n j=1 z j   n j=1 z aj -1 j .
Now we look at the parameters identiability in this case, as it might not hold under the independence assumption.

Proposition 10. The parameters of the ARA 1 model are identiable from the composite likelihood function, that is L (θ|z) = L (θ 0 |z) for all y = (z 1 , . . . , z n ) implies that θ 0 = θ and θ, θ 0 ∈ Θ, as soon as n ≥ 3. In other words, the identiability holds from the likelihood function as soon as observations are conducted three times, and whatever s is.

Proof. Given the expression of the likelihood, we know that the gamma distribution parameters are identiable. Hence, for j = 1, 2, 3, we deduce that

                 αT β = α 0 T β0 b = b 0 (2 -ρ) β = (2 -ρ 0 ) β0 (3 -2ρ) β = (3 -2ρ 0 ) β0
This system is the same one as this in Proposition 8, which was proved to have a unique solution, and so the identiability holds. Now, the log-likelihood is given by

(θ | z) = log(b) n j=1 a j -b n j=1 z j - n j=1 log Γ(a j ) + n j=1 (a j -1) log z j .
The Then, substituting b by its expression above in the log-likelihood leads to

θ | z n =   log   n j=1 a j   -1 -log   n j=1 z j     n j=1 a j - n j=1 log Γ(a j ) + n j=1 (a j -1) log z j
where θ = (α, β, ρ).

We now deal with the extension of this method considering s i.i.d. systems. The log-likelihood has the following expression

(θ | z) = s log b n j=1 a j -s n j=1 log Γ(a j ) + s i=1 log   exp   -b n j=1 z (i) j   n j=1 z (i) j aj -1   6.5. Maximum composite likelihood estimation = s log b n j=1 a j -s n j=1 log Γ(a j ) -b s i=1 n j=1 z (i) j + s i=1 n j=1 (a j -1) log z (i) j .
Once again the parameter b can be expressed with respect to the other ones, leading to

b = s n j=1 a j s i=1 n j=1 z (i) j
and injecting this expression in entails the following expression for the prole composite likelihood function:

θ | z n = log       s n j=1 a j s i=1 n j=1 z i,j       s n j=1 a j -s n j=1 log Γ(a j ) - s n j=1 a j s i=1 n j=1 z i,j s i=1 n j=1 z i,j + s i=1 n j=1 (a j -1) log z i,j
which can be reduced to

θ | z n =   s n j=1 a j     log s -1 + log n j=1 a j -log s i=1 n j=1 z i,j   (6.18) -s n j=1 log Γ(a j ) + s i=1 n j=1
(a j -1) log z i,j .

Composite likelihood based on the increments

As in the previous section, the aim here is to provide the expression of the composite log-likelihood, but dealing with the dierence of the observations

Z j = Z jT --Z (j-1)T -. Let us recall that Z j = W j -V j-1 .
Thus, the p.d.f. of the Z j 's for 1 ≤ j ≤ n is given by f Wj -Vj-1 (z), whose expression is close to that in Equation (6.16), hence we can write

f Zj (z) = b wj +vj-1 Γ(w j )Γ(v j-1 ) exp (-bz) ∞ max(0,-z) (z + x) wj -1 x vj-1-1 exp (-2bx) dx.
Similarly as in the previous case, assuming that the Z j 's are independent allows us to write the joint density of Z, whose expression is given below Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one

f Z (z) = b n j=1 wj +vj-1 n j=1 Γ(w j )Γ(v j-1 ) exp   -b n j=1 z j   n j=1 ∞ max(0,-zj ) (z j + x) wj -1 x vj-1-1 exp (-2bx) dx
and thus the composite log-likelihood has the following expression

(θ | Z) = log(b) n j=1 (w j + v j-1 ) -b n j=1 z j - n j=1 log Γ(w j ) - n j=1 log Γ(v j-1 ) + n j=1 log ∞ max(0,-zj ) (z j + x) wj -1 x vj-1-1 exp (-2bx) dx .
Now we extend the log-likelihood expression to the case where s i.i.d. systems are observed. The joint

density of Z is f Z (z) = s i=1 b n j=1 wj +vj-1 n j=1 Γ(w j )Γ(v j-1 ) exp   -b n j=1 z (i) j   × n j=1 ∞ max(0,-z (i) j ) (z (i) j + x) wj -1 x vj-1-1 exp (-2bx) dx or equivalent f Z (z) = b s n j=1 wj +vj-1 n j=1 Γ(w j ) s Γ(v j-1 ) s exp   -b s i=1 n j=1 z (i) j   s i=1 n j=1 I i,j (z, θ)
where

I i,j (z, θ) = ∞ max(0,-z (i) j ) (z (i) j + x) wj -1 x vj-1-1 exp (-2bx) dx
Finally, the composite log-likelihood is given by

(θ | Z) = s log b n j=1 (w j + v j-1 ) -s n j=1 (log Γ(w j ) + log Γ(v j-1 )) (6.19) -b s i=1 n j=1 z (i) j + s i=1 n j=1 log I i,j (Z, θ).
Once again, the identiability is not studied here for the same reasons as those explained in Section 6.4.2.

Chapter 7. Simulation study α ∈ [0.1, 5] for all methods except the ME method;

β ∈ [0.1, 5] for both ME and CLO methods, the range is reduced to β ∈ [0.1, 2.5] for the HDEI, the HDOI and the CLI methods, due to numerical issues for large beta (see below);

ρ ∈ [0.01, 0.99] for all the methods;

b ∈ [0.1, 5] for the HDEI, the HDOI and the CLI methods.

The numerical issues for large beta are related to the computations of the integrals in the log-likelihood (see Equations (6.15), (6.17) and (6.19)). To be more precise, a large β induces high shape parameters for the involved gamma p.d.f., beyond 3150 when (α, β, ρ) = (1, 5, 0.5) and n = 10, leading to integrals which seem dicult to compute through the standard methods implemented in R. Reducing the maximum possible value of β to 2.5 allows the computations to be possible as the shape parameter does not exceed 28.

Concerning the initialization of the optimization procedure, we have experimented several dierent initial points. In each case, the nal results were almost the same. Then, we have chosen to consider the true values of the parameters as initial points, which allows to slightly reduce the CPU times without aecting the numerical results.

The relative bias (in percentage) and the variances resulting of these estimations are shown in Table 7.1 and Table 7.2 respectively. Also, Figures 7.1 First, from Tables 7.1 and 7.2 as well as Figure 7.5, it can be seen that the CLI method is biased, especially regarding the parameter ρ: the mean bias is beyond 10% associated with a variance less than Regarding the parameters α, β and b, the MC and the QMC methods are better than the MLE, while it is the opposite for the estimation of ρ. We now deal with the cases n = 3, 4, four estimation methods are tested: the MC and QMC based MLE, the CLO and the HDOI method. The other three ones are not considered due to their poor quality. The estimation framework is the same as before except for the number of observed systems, which is reduced to s = 100. Note that even if n < 5 here, the identiability for the HDOI method holds as β = 1 (Corollary 2). The results are given in Table 7.4 Regarding the small size of the data sets (100 trajectories observed at 4 times), the estimations of α and b by the CLO method are rather good. However, this method is not reliable for estimating parameters β and ρ.

In regards to the MC and QMC methods, the rst one is clearly better in terms of bias while both have similar variances. Hence, the MC method is better than the QMC method. Also, unlike the CLO method, the MC method provides good results regarding β and ρ, and the estimation quality of α and b is acceptable. Finally, the HDOI method does not dier for n = 3 and 4 since it considers the rst and third observations only. This method is slightly more ecient for estimating b when n = 3 and α when n = 3 or 4 than the MC method, and both methods have similar variances.

If n = 2, 3, 4, no method stands out clearly from the others.

Let us recall that the MC and QMC approximations deal with the integrals of the function G θ (x, z i,n ) with respect to x over [0, 1] n-1 , for 1 ≤ i ≤ s, whose expression is given by

G θ (x, z i,n ) = n-1 j=1 f Vj C (i) j -m i,j+1 x j + m i,j+1 × n-1 j=1 f Uj z j + C (i) j-1 -m i,j x j-1 -C (i) j -m j+1 x j + m i,j -m i,j+1 × f Wn z i,n + x n-1 C (i) n-1 -m i,n + m i,n
where C j = j i=1 z j and m j = max(0, -z j ). We now study the function G and the integral approximation of this function through the following example.

Example 9. Let us set n = 3 and s = 4. An observations set z 4,3 is generated within the framework of this chapter for the model parameters, leading to the following values: points. However, the mean CPU time required for the retained estimation methods is another point that needs to be enlightened, our CPU times calculations are summarized in Table 7.9.

The mean CPU time for both MC and QMC methods is around 5 to 42 times higher than for the others methods, and this gap is proportional to the number of points for the integral approximations.

Finally, the MC method is preferred to estimate the parameters when associated with a higher number of points, because of its eciency. Also, the use of importance sampling could reduce the CPU time. However, due to their lower CPU times, the MLE method, without using neither MC nor QMC approximations, or the HDOI method can be more suitable to perform further investigation as the estimation of the variance estimators by boostrapping.

All the methods main characteristics are summarized in Table 7.10. Some of the proposed estimation methods can not be used in any case as they are not ecient enough when tested on large samples: the ME, the HDEI and the CLI methods. Also, depending on the number n of observations by trajectory, identiability issues limit the number of possible methods.

As we just stated above, if n < 5, then either the MC method, or the MLE and the HDOI methods are preferred. When n = 10, the HDOI method as well as the CLO method provide good results. In order to determine which one of these two methods is the most ecient, they are tested on various samples of observations in the following section.

Large scale numerical tests

In this part, we deal with the CLO and the HDOI methods, which are tested on various observations sets. The considered model parameters as well as the observations characteristics are partly the same.

To be more precise, the parameters α and b as well as the period of repairs T remain equal to 1, while the other parameters are varying as follows:

Shape function parameter: β = 0.5, 1, 1.5;

Maintenance actions eciency: ρ = 0.2, 0.5, 0. , 5] for the HDOI method. However, the initialization is dierent here: the initial value for the parameters is set at the middle of the search intervals, except in the case ρ = 0.5 where this parameter is initialized at 0.75.

Note that when n = 4 and β = 1.5, the identiability for the HDOI method must be discussed depending on the value of ρ. When ρ = 0.2, 0.5, Corollary 2 states that identiability holds since h 2 (0.2, 1.5) < 1 and h 2 (0.5, 1.5) < 1, but this is not true if ρ = 0.8. However, the numerical study of the identiability, whose results are illustrated in Figure 6.4, led us to conjecture that identiability holds outside a given area, and the point (β, r) = (β, 1 -ρ) = (1.5, 0.2) does not belong to this area. For the sake of clarity, as the RB can be higher than 200% as well as close to zero, instead of representing RB we chose to represent the quantity log(1+RB). Such a scale for the RB is unusual, hence, red dotted lines point out thresholds for which the RB is equal to 5%, 10%, 25% and 50%.

Finally, to be able to compare the results from one parameter to another, as well as the parameters sets, the y-axis range are the same for all plots. However, both the RB and the variance can be higher than the maximum value displayed in the y-axis. This signies that if a curve does not appear on a graph, as the black one in the top right-hand graph of Figure C.5, each point of this curve is higher than this maximum.

As already stated in the previous part, the CLO method cannot be used if n is small because the estimations of β and ρ are biased. This can be seen once again when n = 4 through the gures exposed in Appendix C. In this case, the HDOI method is therefore the most reliable one whatever the value of β and ρ. Moreover, this method provides satisfying results even if only few observations are taken into account.

When n = 16, the value of the parameter ρ does not aect signicantly the methods performance, unlike the parameter β. If β = 0.5, the CLO method appears to provide satisfying results while the HDOI method is biased with variances close to zero, and the value of s has little impact on both the bias and variances. This is clearly due to the number of observations the method deals with, because based on the same observations sets, this method is more ecient considering n = 4. We now try to understand why the performance of the HDOI method becomes poorer when n increases and look at an example.

Example 10. We place ourselves within the following framework: n = 16, T = 1 and θ = (α, β, ρ, b) = (1, 0.5, 0.5, 1). The expectations of the components of Z are given by

E (Z) = ω j (ξ) -ν j-1 (ξ) 1≤j≤16
and plotted in Figure 7.7. This example shows that if the shape function for the underlying degradation process is concave, then the increments of the degradation get smaller and smaller, and vary more and more slowly as j increases.

In Example 10, we can see that a small value for β quickly entails small deterioration increments between PM actions. Let us recall that the maintenance actions put back the system in the situation it was ρT units of times earlier. Hence, because of the slow degradation rate, the maintenance actions have a small impact on the degradation since the system was near the same situation before a repair, whatever its eciency. As a result, this induces a numerical identiability issue, leading to poor quality estimation results. Moreover, the HDOI method only takes into account one out of two increments, which reinforces this problem.

As a result, the rst observations are the most important ones, and considering more observations which contain few additional information about the model parameters adversely aects the estimation quality. However, β = 0.5 leads to a concave shape function and thus the degradation intensity is 7.2. Large scale numerical tests decreasing over time. Therefore, the maintenance actions reduce the deterioration level of the system, but they also induce an increase in the degradation intensity as the system is put back in a state it was earlier, which is inconsistent. As a conclusion, the ARA models are not the most suitable models in this case.

If β = 1, the CLO and the HDOI methods provide similar results for the estimations of α. Regarding b, the rst one is better than the second one, while it is the opposite regarding β and ρ. Finally, when β = 1.5, both methods provide similar results for the estimations of β and ρ, and the CLO method is better in estimating α and b.

Given these characteristics, the choice of the method depends on the context, that is the degradation intensity, the maintenance actions eciency and the number of observations by trajectory. This choice can be done through Table 7.11. Overall, the HDOI method overcomes the CLO method. However, the CLO method always provides good estimations of α and b, whatever the values of β and n. Furthermore, the HDOI method better estimates β and ρ, reducing articially the value of n in some cases. Hence, further work could be to improve the estimation quality by combining both the HDOI and the CLO methods. One approach is to keep the estimations of α and b provided by the CLO method and those of β and ρ from the other one.

In a similar way, plug-in estimations could be tested, that is substituting the values α and b obtained by the CLO method in the log-likelihood of the HDOI method and then estimate β and ρ, as well as the opposite.

Part IV

Conclusion

This work has studied imperfect repair models in the context of gamma deteriorating systems. More specically, estimation methods have been developed for the arithmetic reduction of degradation models of order one and innity (ARD 1 and ARD ∞ ) in Part II, as well as for the arithmetic reduction of age model of order one (ARA 1 ) in Part III.

In Part II, we rst dealt with two classical estimation methods in a fully parametric framework, namely the Moments Estimation and the Maximum Likelihood Estimation (ME and MLE) methods.

The parameters identiability was studied, estimators expressions were provided and nally the performance of the methods was illustrated. After that, arising from the study of the MLE method, an estimator for the maintenance actions eciency (parameter ρ) was proposed in a semiparametric framework. In the case of one single trajectory and when the shape function of the underlying gamma process is concave, this estimator was proved to be strongly consistent as the number of repairs tends towards innity, with a surprisingly high convergence rate, at least exponential for some particular cases. This work was next extended to the case where s independent and identical systems are observed. A similar semiparametric estimator was proposed for ρ and the strong consistency was proved to hold as s tends towards innity, no matter the xed number of repairs and out of any technical condition requirement. The convergence rate was also studied, which was shown to depend on the shape function of the gamma process and on the maintenance period, leading to a speed that can be either slower or faster than √ s, according to the case.

In Part III, we proposed several estimations methods for an ARA 1 model, which are based on either the observations or on the increments, leading to six dierent estimation methods. Depending on the method, either expressions of the estimators or the log-likelihood are provided, and the identiability of the model parameters was studied regarding four out of the six methods. In order to study and compare the performance of all these methods, numerical investigations based on simulated data have been conducted at a large scale. Two of them has appeared to be more ecient than the other ones: the Half Data based on the Odd Indexes (HDOI) and the Composite Likelihood based on the Observations (CLO) methods. Overall, however, none of these two methods stands out clearly from the other.

Thus the choice of the method must be done with respect to the context. To be more specic, this choice depends on the number of repairs, the repairs eciency and the degradation intensity evolution over time.

Several points of interest would be interesting to study in complement to this thesis. There exist many other models that extend the ARD models studied here, such as for instance the ARD m model for which the basic idea is that a maintenance action removes a proportion of the degradation accumulated by the system from the last m maintenance actions. An idea could be to generalize the semiparametric estimate of the maintenance eciency to such a model, which seems to be possible as this method relies on the non negativity of the gamma process. Moreover, the adaptation of this estimation procedure to another monotonous Lévy process than the gamma process would be interesting to study. Finally, based on a short numerical study, we have seen that the condition which ensures an exponential convergence rate was sucient but not necessary. Further investigations could be done in order to see whether it could be possible to rene the mathematical conditions under which the dierent convergence rates are obtained.

Regarding the work done in Part III, several estimation procedures could be improved. As already stated, it could be interesting to see whether the estimation quality is better by combining the HDOI and the CLO methods. Moreover, we explained that the Monte Carlo and Quasi Monte Carlo integral approximations required in the MLE method could be improved by using importance sampling. Further work should be done to address these possible improvements.

Beyond that, except in Chapter 4, asymptotic properties of the estimators were not investigated in this document and further work is required for a better understanding of their behaviour, beyond the numerical investigation performed in this thesis. One might study these properties as an extension of the present work. Besides all of this, further work could be the development of estimation procedures in the context of an ARA ∞ model, for which a maintenance action removes a proportion of the age accumulated by the system since it was put into operation. More research might be done in this context, mimicking the work done regarding the ARA 1 model as a rst step.

Finally, another point could be to consider an observation scheme which is decoupled from the maintenance schedule. It is not the case here as the degradation level is measured right before each repair.

As an example, observations could still be conducted periodically while the system could be maintained at random times according to a Poisson process. However, preliminary investigations were done in this framework, and the degradation level of the maintained system seems dicult to write down with a simple expression, but that could remain exploitable for developing estimation procedures. Hence, this seems a challenging (but interesting) subject for further research. et

Y 2T = Y 2T -+ ρ X (2)
2T --X

(2)

T = (1 -ρ) X (2)
2T -X

(2) T + X

T -X

(1) 0

.

Ainsi, pour tout t dans [ nT , (n + 1)T [, avec n dans N * , on a :

Y nT = (1 -ρ) n j=1 X (j)
jT -X

Y t = X (1) t et Y T = (1 -ρ) X (1)
T -X

(1) 0

= (1 -ρ)X (1)
T , et pour tout t dans [ T , 2T [

Y t = Y T + X (2) t -X (2) T
.

Lorsque la seconde maintenance survient, une proportion ρ de la dégradation accumulée depuis l'instant t = 0 est retirée, ainsi

Y 2T = (1 -ρ) Y T + X (2)
2T -X

(2) T

, qui peut être écrit comme

Y 2T = (1 -ρ) 2 X (1)
T -X

(1) 0

+ (1 -ρ) X (2)
2T -X

(2) T .

  Figure D.1). Many models have been envisioned in the literature for their Time Failure intensity Failure intensity Minimal maintenance action Imperfect maintenance action Perfect maintenance action

Figure 1

 1 Figure 1

Figure 2 :

 2 Figure 2: Simulated trajectories of a gamma process with shape function a(t) = t and scale parameter b = 1.

1. 2 .Figure 1 . 1 :

 211 Figure 1.1: Expectations of X t and Y t , where the shape function is such that a : t → αt β with α = 2 and β = 3/2, and with b = 1 and ρ = 1/2.

Figure 1 .Figure 1 . 2 :

 112 Figure 1.2: Variances of X t and Y t , where the shape function is such that a : t → αt β with α = 2 and β = 3/2, and with b = 1 and ρ = 1/2.

,Figure 1 . 3 :

 13 Figure 1.3: Expectations of X t and Y t , where the shape function is such that a : t → αt β with α = 2 and β = 3/2, and with b = 1 and ρ = 1/2.

Figure 1 . 4 :

 14 Figure 1.4: Variances of X t and Y t , where the shape function is such that a : t → αt β with α = 2 and β = 3/2, and with b = and ρ = 1/2.

Chapter 2 .

 2 Parametric inference for the Arithmetic Reduction of Degradation models Proof. The probability multiplication rule provides f Y (y) = f Y1 (y 1 ) n j=2

  ) and(2.18) provide an estimate θ of the model parameters.

Chapter 2 .

 2 Parametric inference for the Arithmetic Reduction of Degradation models The aim of the next chapter is to investigate numerically the estimation methods performances for each model, starting with the ARD 1 model. Reduction of Degradation model of order one We place ourselves in the following framework for the model parameters and the observations characteristics: Shape function parameter: α = 1 and β ∈ {1.2, 1.6}; Scale parameter: b = 1; Maintenance actions eciency: ρ ∈ {0.2, 0.5, 0.8}; Period of repairs: T = 1; Observations times: {jT -; 1 ≤ j ≤ n} with n ∈ {5, 20}; Number of observed i.i.d. systems: s ∈ {50, 200}.

  Figure 4.1 for a(t) = t 1.5 , b = 1, T = 1 and ρ = 0.5, together with the corresponding trajectories of the (X (j) t ) t≥0 , j = 1, 2, . . .

Figure 4 . 1 :

 41 Figure 4.1: An example of simulated trajectories of X (j) t t≥0 , j = 1, . . . , 5 and of (Y t ) t≥0 (intrinsic and overall degradation levels, respectively) with parameters a(t) = t 1.5 , b = 1, T = 1 and ρ = 0.5.

Proposition 6 .

 6 (Identiability) Let Y = Y 1 , . . . , Y n and Ỹ = Ỹ1 , . . . , Ỹn be two random vectors based on the ARD 1 repair model with parameters (a(•), b, ρ) and ã(•), b, ρ , respectively (and the same period T ). Assume that there are at least two observations (n ≥ 2). Then if Y and Ỹ are identically distributed (denoted by Y D = Ỹ ), necessarily ρ = ρ. Proof. Assume that Y D = Ỹ and n ≥ 2. Then Y 1 and Ỹ1 are identically distributed, with Y 1 = U 1 ∼ Γ (a 1 , b) and Ỹ1 = Ũ1 ∼ Γ ã1 , b . This implies that a 1 = ã1 and b = b. Also, Y 2 = (1 -ρ)U 1 + U 2 and Ỹ2 = (1 -ρ) Ũ1 + Ũ2 must share the same distribution, and hence the same Laplace transform.

Figure 4 . 2 :

 42 Figure 4.2: Representation of the sets D j , j = 1, . . . , n (black horizontal segments) for ρ 0 = 0.5 (vertical blue line) and b = 1, with a(t) = √ t (concave function), n = 30 for the left plot, and a(t) = t 1.5 (convex function), n = 10 6 for the right one, where the lower bounds M j 's of the D j 's are highlighted by blue crosses, Example 2.

For a better

  understanding of what the D n 's are, we now look at an example, based on simulated data, where we consider the successive D 1 , . . . , D n (where the D j 's, j = 1, . . . , n, are dened in a similar way as D n ).

(4. 8 )Lemma 4 .

 84 Let us denote by I x (α 1 , α 2 ) the cumulative density function of the beta distribution B (α 1 , α 2 )

(4. 16 )Chapter 4 .

 164 Semiparametric estimate of the maintenance actions eciency with u n = C ε a 2(n+1) 4n

4. 4

 4 Extension to the case where several systems are observed 4.4.1 Extended semiparametric estimator In this section, s identical and independent systems are considered. They share the same intrinsic deterioration and ARD 1 repair model with parameter (a(•), b, ρ 0 ) and they are all observed at times T -,

Thus, M s, 2

 2 tends towards ρ 0 almost surely, which proves the result. Theorem 4. Let n ≥ 2. The almost sure convergence rate of the estimator M s,n (with respect to s) is at least s -k , for all positive real number k such that

1≤i≤s ρ 0

 0 -M (i) n > εs = P ρ 0 -M (1) n > εs s ≤ P min 2≤p≤n Up U p-1 > εs s , which leads to the result from Corollary 1.Chapter 4. Semiparametric estimate of the maintenance actions eciency convergence of the series with generic term u s . If 1 -ka 2p ≤ 0, then u s converges towards 1 or exp(-C p ) (if ka 2p = 1), and the series is divergent. If 1 -ka 2p > 0, then lim s→+∞ s 2 u s = 0 and u s = s→+∞ o 1

1 ; 5 ; 1 ;

 151 Shape function: a(t) = log(1 + t) (ECR) and a(t) = √ t (S-ECR); Scale parameter: b = Maintenance eciency parameter: ρ 0 = 0.Period for repairs: T = Observation times: {nT -; 1 ≤ n ≤ 250} because the system is maintained 250 times. Thus for a single maintained system simulated over the time interval [0,250], that is a degradation trajectory over the time interval [0,250], M j is computed for each observation time (right before the repair time) providing a realization (m 1 , . . . , m 250 ) of the random vector (M 1 , . . . , M 250 ). We generate 250 000 i.i.d trajectories, which leads to 250 000 i.i.d. realizations m (i) 1 , . . . , m(i)250 ; 1 ≤ i ≤ 250 000 of (M 1 , . . . , M 250 ). In other words, we have 250 000 estimations of ρ 0 at each observation time T -, 2T -, . . . , 250T -. Then the Empirical Bias (EB) is computed at each observation time nT -as follows

Figure 4 . 3 :

 43 Figure 4.3: Plots of the common logarithm of the Empirical Bias versus n, for the ECR case on the left and the S-ECR case on the right.

Figure 4 . 4 :

 44 Figure 4.4: Plots of the PNEE as a function of n, for the ECR case on the left and the S-ECR case on the right. Blue vertical lines point to the smallest repair numbers n for which PNEE(n) is larger or equal than 95%, 99.5%, 99.99% and 100% from left to right.

Figure 4 . 5 :

 45 Figure 4.5: Partial plots of the common logarithm of the Empirical Bias versus n, for the ECR case on the left and the S-ECR case on the right. The blue lines correspond to the rst degree polynomial regression for the ECR case and to the second degree polynomial regression for the S-ECR case.

Figure 4 . 6 :

 46 Figure 4.6: Plots of the common logarithm of the Empirical Bias versus n, with a 2n = n -1 log n

Figure 5 . 1 :

 51 Figure 5.1: Expectations (equal to variances) of X t and Z t , where the shape function is such that a : t → αt β with α = and β = 3/2, and with b = 1 and ρ = 1/2.

  with D(θ) the empirical version of D(θ , θ 0 ) (see Equation (2.1)). Let us set n ≥ 3 and d = 2, and recall the used parametrization for the ME method: θ = (µ, η, β, ρ) with µ = α/b and η = α/b 2 , that is α = µ 2 /η and b = µ/η). Thus, we have

γ 0 ≤ 1

 1 and rβ r0β0 ∈ (w (1) , w (2)) with w (1) < w (2); γ 0 ≥ 2 and rβ r0β0 ∈ (w (2) , w (1)) with w (2) < w (1); γ 0 ∈ (1, 2) and rβ r0β0 ∈ (w (γ 0 ) , max (w (1) , w (2))).

Corollary 4 .

 4 Suppose Conditions (*) to hold. If the function u 2 is non-monotonic, then γ 0 < 2 and rβ r0β0 ∈ (w (1) , w (2)), with w (1) < w (2). Proof. From the second point of Corollary 3, if u 2 is non-monotonic and γ 0 ≥ 2, then rβ r0β0 ∈ (w (2) , w (1)).

Figure 6 . 1 :

 61 Figure 6.1: Both graphs highlight the possible areas for (r, β) fullling Conditions (*) given a xed (r 0 , β 0 ), based on Corollary 4. The blue hatched area corresponds to the (r, β) such that γ 0 < 2 and rβ r0β0 > w(1), while the grey hatched area corresponds to the (r, β) such that γ 0 < 2 and rβ r0β0 < w(2). The red crosses represent the values of (r 0 , β 0 ), with (r 0 , β 0 ) = (0.01, 5) on the left-hand plot and (r 0 , β 0 ) = (0.05, 4) on the right-hand one.

  rβ r0β0 ≈ 1.03 ∈ w(1), w(2) where w(1) ≈ 1.02 and w(2) ≈ 1.05.

Figure 6 . 2 :

 62 Figure 6.2: The surface represents the common logarithm of g (r0,β0) (r, β) depending on (r, β), with (r 0 , β 0 ) ≈ (0.07, 4.14). Both blue lines highlight the two points (r 0 , β 0 ) and (r, β) = (0.1, 3) for which g (r0,β0) (r, β) ≈ 0.

Chapter 6 .Figure 6 . 3 :

 663 Figure 6.3: The possible area for (r, β) fullling Conditions (*) given (r 0 , β 0 ) ≈ (0.07, 4.14), based on Corollary 4, is plotted in blue. Between the curves, we have γ 0 < 2 and rβ r0β0 ∈ (w (1) , w (2)), with w (1) < w (2). The red crosses represent (r, β) and (r 0 , β 0 ).

  Figure 6.4: The black points are the values for (r 0 , β 0 ) for which there exists (r, β) such that

Figure 6 . 5 :

 65 Figure 6.5: Both areas match the points (r, β) for which either the results of Corollary 4 (blue area) or the results of Corollary 5 (grey area) are satised, with β 0 > β, r > r 0 and (r 0 , β 0 ) ≈ (0.07, 4.14).

  parameter b can be written as a function of the other ones by solving the equation ∂ b (θ | z) = 0,

  to 7.5 represent the pairs plots of the estimations for each of the ve methods.The classical MLE is not treated here because it requires the use of MC and QMC approximations, and the large dimension (n -1 = 9) induces too high computing times.

Figure 7 . 2 :

 72 Figure 7.2: Pairs plots of HDEI estimates. The red lines indicate the true value of each parameter.

7. 1 .

 1 Methods selection

Figure 7 . 3 :

 73 Figure 7.3: Pairs plots HDOI estimates. The red lines indicate the true value of each parameter.

Figure 7 . 4 :

 74 Figure 7.4: Pairs plots of CLO estimates. The red lines indicate the true value of each parameter.

Figure 7 . 5 :

 75 Figure 7.5: Pairs plots of CLI estimates. The red lines indicate the true value of each parameter.

7. 1 .Figure 7 . 6 :

 176 Figure 7.6: Levels of the function G θ0 (x 1 , x 2 , z i,3 ) depending on x 1 , x 2 and for i = 1, 2, 3, 4 from left to right, top to bottom.

8 ;

 8 Observations times: {jT -; 1 ≤ j ≤ n} with n = 4, 16; Number of observed i.i.d. systems: s = 5, 25, 100, 400 and the method is used 2000 times for each combination of those parameters. The estimations are computed as before, that is a maximization by the gradient method, searching (α, β, ρ) over [0.1, 5] 2 × [0.01, 0.99] for the CLO method, and (α, β, ρ, b) over [0.1, 5] × [0.1, 2.5] × [0.01, 0.99] × [0.1

Figure 7 . 7 :

 77 Figure 7.7: Expectations of the degradation increments with respect to j. The red dots matches the information taken into account in the HDOI method while the grey ones are not.

2 Figure A. 1 : 2 Figure A. 2 : 5 Figure A. 3 : 5 Figure A. 4 : 8 Figure A. 5 : 8 Figure A. 6 : 2 Figure A. 7 : 2 Figure A. 8 : 5 Figure A. 9 : 5 Figure A. 10 : 8 Figure A. 11 : 8 Figure A. 12 : 2 Figure B. 1 : 2 Figure B. 2 : 5 Figure B. 3 : 5 Figure B. 4 : 8 Figure B. 5 : 8 Figure B. 6 : 2 Figure B. 7 : 2 Figure B. 8 : 5 Figure B. 9 : 5 Figure B. 10 : 8 Figure B. 11 : 8 Figure B. 12 :Figure C. 1 :Figure C. 3 :Figure C. 5 :Figure C. 7 :Figure C. 9 :Figure C. 11 :Figure C. 13 :Figure C. 15 :Figure C. 17 :

 2122535485862728595108118122122535485862728595108118121357911131517 Figure A.1: Plots of the RB evolution of the parameters estimations depending on s, when (β, ρ) = (1.2, 0.2).

  nT suivent des lois gamma Γ (a(nT ) , b/(1 -ρ)) et Γ (a(t) -a(nT ) , b) respectivement.Concernant le modèle ARD ∞ , une action de maintenance à l'instant jT réduit le niveau de dégradation à cet instant de ρ%, soit la dégradation du système accumulée depuis l'instant initial t = 0. Ce modèle et le précédent se comportent de manière identique sur l'intervalle de temps [0 , 2T [. Ainsi, en conservant les mêmes notations que pour le modèle précédent, on peut directement écrire que pour tout t dans [0 , T [ on a

  Chapter 4. Semiparametric estimate of the maintenance actions eciency Now, assume Points 3-5 to be true for some j ∈ {1, • • • , n -1} and, to begin with, let us note that Equation (4.4) implies that

). Hence Points 3-5 are true for j = 1.

Table 4 .

 4 1: Linear regressions summary for the ECR case (value / p-value) Intercept coecient First degree coecient Second degree coecient 1 st degree 0.10 / 6.12 × 10 -4 -0.36 / < 2 × 10 -16 2 nd degree 0.19 / 1.70 × 10 -5 -0.38 / < 2 × 10 -16 6.9 × 10 -4 / 4 × 10 -3 st degree -0.95 / 4.79 × 10 -13 -0.19 / < 2 × 10 -16 2 nd degree -0.21 / 4.90 × 10 -6 -0.27 / < 2 × 10 -16 1.76 × 10 -4 / < 2 × 10 -16

	(Residual error)	Minimum	1 st quartile	Median	3 rd quartile	Maximum
	1 st degree	-0.09	-0.05	-0.01	0.04	0.12
	2 nd degree	-0.14	-0.03	3 × 10 -3	0.04	0.10
		Table 4.2: Linear regressions summary for the S-ECR case
	(value / p-value)	Intercept coecient	First degree coecient Second degree coe-
					cient
	1 (Residual error)	Minimum	1 st quartile	Median	3 rd quartile	Maximum
	1 st degree	-0.31	-0.27	-0.14	0.21	0.84
	2 nd degree	-0.15	-0.06	-5 × 10 -3	0.06	0.22

Table 7 .

 7

			1: Relative bias estimation
	Methods / Relative bias (%)	α	β	ρ	b
	ME	0.14%	1.18%	2.43%	0.01%
	HDEI	1.32%	0.11%	0.06%	0.63%
	HDOI	0.11%	0.002%	0.17%	0.1%
	CLO	0.08%	0.29%	0.02%	0.01%
	CLI	2.79%	1.91%	10.24%	2.27%

Table 7 .

 7

	2: Variances estimation

Table 7 .

 7 3: Relative bias and variances for the MLE, MC and QMC methods when n = 2

	Methods / Relative bias (%)	α	β	ρ	b
	MLE	2.66%	1.12%	0.62%	3.05%
	MC	0.08%	0.52%	7.01%	0.03%
	QMC	0.45%	0.43%	3.52%	0.77%
	Methods / Variances	α	β	ρ	b
	MLE	9.7 × 10 -3	2.2 × 10 -2	4.7 × 10 -3	1.2 × 10 -2
	MC	1.7 × 10 -2	3.7 × 10 -2	1.6 × 10 -2	2.0 × 10 -2
	QMC	1.5 × 10 -2	3.9 × 10 -2	1.5 × 10 -2	2.2 × 10 -2

Table 7 . 4 :

 74 to Table7.7. Relative bias and variances for the CLO method when n = 3, 4

	Relative bias (%)	α	β	ρ	b
	n = 3	0.43%	42.56%	3.44%	0.39%
	n = 4	0.29%	25.44%	6.36%	0.13%
	Variances	α	β	ρ	b
	n = 3	5.7 × 10 -3	0.62	0.11	1.1 × 10 -3
	n = 4	5.3 × 10 -3	0.39	0.08	7.0 × 10 -4

Table 7 .

 7 5: Relative bias and variances for the MC method when n = 3, 4

	Relative bias (%)	α	β	ρ	b
	n = 3	2.98%	0.12%	1.46%	2.80%
	n = 4	6.65%	0.28%	0.43%	4.34%
	Variances	α	β	ρ	b
	n = 3	1.7 × 10 -2	2.6 × 10 -2	8.6 × 10 -3	1.8 × 10 -2
	n = 4	8.8 × 10 -3	9.4 × 10 -3	3.3 × 10 -3	1.3 × 10 -2

Table 7 . 6 :

 76 Relative bias and variances for the QMC method when n = 3, 4

	Relative bias (%)	α	β	ρ	b
	n = 3	6.14%	4.93%	3.19%	7.27%
	n = 4	5.01%	1.92%	3.17%	4.57%
	Variances	α	β	ρ	b
	n = 3	1.5 × 10 -2	2.8 × 10 -2	7.1 × 10 -3	1.5 × 10 -2
	n = 4	1.1 × 10 -2	1.1 × 10 -2	5.5 × 10 -3	1.0 × 10 -2

Table 7 .

 7 7: Relative bias and variances for the HDOI method when n = 3, 4

		α	β	ρ	b
	Relative bias (%)	2.22%	2.19%	1.73%	3.54%
	Variances	1.3 × 10 -2	1.7 × 10 -2	8.5 × 10 -3	1.8 × 10 -2

Table 7 .

 7 9: Mean CPU times for one single estimation when n = 2, 3, 4, based on 500 repetitions

		MC and QMC	MLE	HDOI
	n = 2	10min	2min	NA
	n = 3	16min	NA	28s
	n = 4	20min	NA	28s

  For each (β, ρ) in {0.5, 1, 1.5} × {0.2, 0.5, 0.8}, the Relative Bias (RB) and the variances of the estimations are summarized in Appendix C. The gures deal either with the RB or with the variance, and are composed of four graphs: one for each parameter of the model. These graphs show four curves each, representing either the relative bias or the variance evolution with s for both the HDOI and the CLO methods when n = 4 and n = 16.

	Chapter 7. Simulation study			7.2. Large scale numerical tests
			Table 7.10: Characteristics summary for the methods selection
			Numerical ap-
	Methods	Identiability	proximation of	Estimations	Pros		Cons
			integrals
					No	integral
	ME	n ≥ 3	None	Unbiased and high dis-persion	approxima-tions and quick	Lack of accuracy
					computation
	MLE	n ≥ 2	Monte-Carlo and Quasi Monte-Carlo approximations (integrals of di-mension n -1), or numerical ap-proximations of one-dimensional integrals when n = 2	Unbiased only if the number of points for the integral approx-imation is sucient, medium dispersion no matter the number of points. The MC based estimation is better	Identiability holds as soon as n ≥ 2	The number of points for the MC and QMC ap-proximations does not af-fect the estimators' vari-ances, but there a bias appears if it is too low. Unusable when n is too large because the number of points must be larger as well, which leads to non manageable comput-ing time.
			Numerical ap-	Unbiased and medium
	HDEI	Not studied	proximations of one-dimensional	dispersion. Some es-timates are very far	Accurate overall	It is possible for an estim-ation to be an outlier
			integrals	from the true value
	HDOI	n ≥ 5 (par-tially when n = 3, 4)	Numerical ap-proximations of one-dimensional integrals	Unbiased and low dis-persion	Accurate reliable method and	Identiability holds over the entire parameters set only if n ≥ 5
	CLO	n ≥ 3	None	Unbiased and medium to low dispersion	No approxima-integral tions and quick computation, accurate and reliable method	The parameter β can-not be estimated properly (biased and high disper-sion) when n is small, that is n = 3, 4.
	CLI	Not studied	Numerical ap-proximations of one-dimensional integrals	Highly biased (from 2% for α to 10% for ρ) and low dispersion	None		A high bias associated with a low dispersion make this method unus-able

Table 7 .

 7 11: Summary for the choice of the estimation method

		Decreasing degradation in-	Constant degradation intens-	Increasing degradation in-
		tensity (β < 1)	ity (β = 1)	tensity (β > 1)
	n = 4		HDOI
	n = 16	CLO if ρ is small, otherwise HDOI reducing n to 4	HDOI	CLO

Due to the introduction of a published paper within this chapter, some notations dier from those used in the remainder of the document.

ε a 2(n+1) 2 1+a2n+1+a 2(n+1) n

Contents

Part III Imperfect repairs based on reduction of age Chapter 6. Parametric inference for the Arithmetic Reduction of Age model of order one

j is the observed increments of the degradation of the ith maintained system between the instants right before the j -1th and the jth maintenance for 1 ≤ j ≤ n and 1 ≤ i ≤ s, which is a realisation of the r.v. Z j ;

.

We provide another representation of such transformed observations in the lemma below.

Lemma 7. Let Z j be the degradation increment between the instants right before the j -1th and the jth maintenance. For all 1 ≤ j ≤ n, the r.v. Z j can be expressed with respect to independent and gamma distributed random variables as follows:

with the convention V 0 = 0 and where U 1 , . . . , U n ,V 1 , . . . , V n are independent;

and the shape functions µ j and ν j are such that

Moreover, because of the independence of U j and V j , we also can write Z j as

where W j and V j-1 are independent, and W j is gamma distributed Γ (ω j (ξ) , b) with

Proof. Let U j be the increment of the underlying degradation process over the real time interval ((j -1) T, (j -ρ) T ].

Therefore, we have

which can be written as Chapter 7 Simulation study

Methods selection

The aim of this chapter is to investigate numerically the estimators quality for each method. For sake of simplicity, the estimation methods are abbreviated as follows:

Moments estimation: ME;

Maximum likelihood estimation: MLE;

Maximum likelihood estimation using Monte Carlo approximation: MC;

Maximum likelihood estimation using Quasi Monte Carlo approximation: QMC; Maximum likelihood estimation using half of the data which have an odd index: HDOI;

Maximum likelihood estimation using half of the data which have an even index: HDEI;

Maximum likelihood estimation using the composite likelihood based on the observations: CLO;

Maximum likelihood estimation using the composite likelihood based on the increments: CLI;

The estimation methods are rst tested on large samples in order to select the most ecient ones. To be able to compare results, we place ourselves within the same framework for each method. The model parameters and the observations characteristics are the following:

Shape function parameter: ξ = (α, β) = (1, 1);

Scale parameter: b = 1;

Maintenance actions eciency: ρ = 0.5;

Period of repairs: T = 1;

Observations times: {jT -; 1 ≤ j ≤ n} with n = 10; Number of observed i.i.d. systems: s = 1000.

We generate 500 observations sets, and for each method except for MC and QMC an estimation of θ = (α, β, ρ, b) is computed for each set. These estimations are based on the minimization or maximization of the quantities provided in Equations (6.4), (6.15), (6.17), (6.18) and (6. [START_REF] Giorgio | A new state-dependent degradation process and related model misidentication problems[END_REF]), based on a gradient method.

The range for the parameters in the optimization procedure is the following:

7.1. Methods selection 5.5 × 10 -5 . This disqualies the CLI method. Among the remaining estimation methods, the ME method is globally less ecient than the other three:

regarding the parameters β and ρ, the bias and the variances are higher than for any other method. This can also be seen by comparing We observe the same characteristics for the HDEI regarding α and b, that is a higher bias and variance than for the HDOI and CLO methods. Moreover, extreme values can be seen in Figure 7.2. This is a lack of robustness of the HDEI method compared to alternative methods since for some data sets estimations can be very far from the true parameters values. In conclusion, the HDOI and the CLO methods are the most reliable with respect to the three other ones (ME, HDEI and CLI). Hence, in the following, these methods are investigated and compared in order to determine which one is the most ecient. Remember that the identiability holds for n ≥ 5

for the HDOI method while it holds as soon as n ≥ 3 for the CLO method. Moreover, if n = 2, the identiability holds for the MLE only. Hence, before going on further with the comparison between the HDOI and the CLO methods, the cases n = 2, 3, 4 are studied.

As already mentioned, the MC and QMC methods require large computing times. However, in order to increase the estimations accuracy, the number of points has to be large regarding n and this aects the computing time. Hence, the number of systems s and the number of points for the integral approximations are set to 100 and 2 14 = 16, 384 respectively.

Chapter 7. Simulation study

where the component z i,j is the degradation level of the i-th system at time jT -. For each

The level sets of this function are plotted in Figure 7.6 for i = 1, 2, 3, 4 from left to right, top to bottom. It can be seen that there is a region where the function G is zero, which highly depends on the observations and can reach more than half of the integration region. Regarding the remaining area where G is non-zero, the shape of the function is similar in each plot, that is:

slowly increasing from the light grey area, where it is minimal, to the darker grey area;

quickly increasing towards innity at the edge (red area).

Within this framework, the shape parameter of the U j 's and V j 's is 0.5. Hence, their p.d.f., which are considered in the expression of G, are equal to the function

Regarding the r.v. U 1 , we have

). Hence, when x 1 tends towards 1 then y tends towards zero. In this case f U1 tends to ∞, this explains why G tends towards ∞ at this edge, which also can be seen in Figure 7.6. Similar results are obtained regarding the other distributions, which leads to the observed growth over almost all the edges. Note that this particularity is not an issue regarding the p.d.f. of W n since its shape parameter is 1 here. Now, we are interested in the importance of the edge in the integral approximation. With this aim, the integrals are computed with roughly 5 × 10 5 points for the MC. The proportion of points belonging to the red area as well as the proportion of the integral arising from the red area are summarized in Table 7.8. The weight of the edges is high because, for each computation, 6% to 19% of the integral arise from the red area while this area represents less than 1% of the integral domain.

The previous example leads us to the conclusion that the MC and QMC based methods could be improved by using importance sampling. Another alternative could be to consider a larger number of Appendix A ARD 1 model, results of the simulation study As already stated in Section 3, the results are exposed as follows.

For each possible combination of (β, ρ), the RB and the variances of the estimations are summarized in the following. The gures deal either with the RB or with the variance, and are composed of four graphs: one for each parameter of the model. These graphs show four curves each, representing either the relative bias or the variance evolution with s for both the ME and the MLE methods when n = 5

and n = 20.

To be able to compare the results from one parameter to another, as well as the parameters sets, the y-axis range is the same for all plots. However, both the RB and the variance can be higher than the maximum value displayed in the y-axis. This signies that if a curve does not appear on a graph, as the black one in the top right-hand graph as well as the bottom left-hand one of For each possible combination of (β, ρ), the RB and the variances of the estimations are summarized in the following. The gures deal either with the RB or with the variance, and are composed of four graphs: one for each parameter of the model. These graphs show four curves each, representing either the relative bias or the variance evolution with s for both the ME and the MLE methods when n = 5

and n = 20.

To be able to compare the results from one parameter to another, as well as the parameters sets, the y-axis range is the same for all plots. However, both the RB and the variance can be higher than the maximum value displayed in the y-axis. This signies that if a curve does not appear on a graph, as the black one in the top right-hand graph as well as the bottom left-hand one of For each (β, ρ) in {0.5, 1, 1.5} × {0.2, 0.5, 0.8}, the Relative Bias (RB) and the variances of the estimations are summarized in the following. The gures deal either with the RB or with the variance, and are composed of four graphs: one for each parameter of the model. These graphs show four curves each, representing either the relative bias or the variance evolution with s for both the HDOI and the CLO methods when n = 4 and n = 16.

For the sake of clarity, as the RB can be higher than 200% as well as close to zero, instead of representing RB we chose to represent the quantity log(1+RB). Such a scale for the RB is unusual, hence, red dotted lines point out thresholds for which the RB is equal to 5%, 10%, 25% and 50%.

Finally, to be able to compare the results from one parameter to another, as well as the parameters sets, the y-axis range are the same for all plots. However, both the RB and the variance can be higher than the maximum value displayed in the y-axis. This signies that if a curve does not appear on a graph, as the black one in the top right-hand graph of Figure C.5, each point of this curve is higher than this maximum. . Dans la littérature, un grand nombre de modèles ont été envisagés pour les modéliser, tels que, par exemple, les modèles d'âge virtuel introduits par Kijima [START_REF] Kijima | Some results for repairable systems with general repair[END_REF], traités de façon plus approfondie dans [START_REF] De Toledo | ARA and ARI imperfect repair models: Estimation, goodness-of-t and reliability prediction[END_REF]16], et davantage dans [START_REF] Brenière | Virtual age models with time-dependent covariates: A framework for simulation, parametric inference and quality of estimation[END_REF] où les auteurs ajoutent des covariables au modèle d'âge virtuel. D'autres possibilités sont les processus géométriques [START_REF] Lam | The geometric process and its applications[END_REF] (plus largement étudié dans [START_REF] Bordes | Extended geometric processes: semiparametric estimation and application to reliability[END_REF], et plus récemment dans [START_REF] Dauxois | Semiparametric inference for an extended geometric failure rate reduction model[END_REF]) ou, comme évoqués plus haut, les modèles basés sur la réduction de l'intensité de défaillance [START_REF] De Toledo | ARA and ARI imperfect repair models: Estimation, goodness-of-t and reliability prediction[END_REF]16]. Voir, par exemple, [START_REF] Doyen | On geometric reduction of age or intensity models for imperfect maintenance[END_REF] pour un compte rendu récent et des extensions pour de tels modèles. Voir aussi [START_REF] Pham | Imperfect maintenance[END_REF] De nos jours, l'évolution de ce nombre réel au cours du temps est communément modélisée par un processus stochastique, souvent considéré comme ayant une tendance positive. Les modèles classiques incluent le processus inverse Gaussien [START_REF] Ye | The Inverse Gaussian Process as a Degradation Model[END_REF] ou les processus de Wiener (avec tendance) [START_REF] Hu | Optimum step-stress accelerated degradation test for wiener degradation process under constraints[END_REF][START_REF] Liu | A condition-based maintenance policy for degrading systems with age-and state-dependent operating cost[END_REF][START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF]. Dernièrement, le processus de Wiener transformé a également été introduit par [START_REF] Giorgio | A new age-and state-dependent degradation process with possibly negative increments[END_REF], et davantage étudié dans [START_REF] Giorgio | Bayesian estimation and prediction for the transformed Wiener degradation process[END_REF], où les incréments de la dégradation peuvent être négatifs. Tous ces processus stochastiques sont relativement commun dans des domaines variés autre que celui de la abilité, tel que la nance, les assurances ou l'épidémiologie. Cette thèse est dédiée au processus gamma, largement répandus depuis qu'ils ont été introduits dans le domaine de la abilité par Abdel-Hameed [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF] et Çinlar [START_REF] Çinlar | Stochastic process for extrapolating concrete creep[END_REF]. Ce processus étant monotone, il est parfaitement adapté pour la modélisation de l'évolution de dégradation croissante.

Avant de dénir le processus gamma, la loi gamma est dénie, permettant d'introduire les notations et la paramétrisation utilisées dans la suite du document.

Une variable aléatoire X suit une loi gamma de paramètre de forme a > 0 et de paramètre d'échelle b > 0 (noté X ∼ Γ(a, b)), si sa fonction de densité est telle que :

par rapport à la mesure de Lebesgue. L'espérance et la variance d'une telle variable sont données par

On rappelle maintenant quelques propriétés très connues de la loi gamma : Soient X 1 et X 2 deux variables aléatoires indépendantes suivant des lois gamma Γ(a 1 , b) et Γ(a 2 , b) respectivement, avec a 1 , a 2 , b > 0. Ainsi, pour tout c strictement positif, les variables aléatoires cX 1 et

On pose maintenant la fonction a(•) : R + → R + croissante et continue, telle que a(0) = 0, et on pose b > 0. Soit (X t ) t≥0 un processus stochastique continue à droite avec limites à gauche. Le processus (X t ) t≥0 est un processus gamma non homogène ayant pour fonction de forme a(•) et pour paramètre d'échelle b, si X 0 = 0 presque sûrement; les incréments du processus sont indépendants; les incréments suivent des lois gamma, autrement dit pour tout 0 ≤ s < t, on a X t -X s ∼ Γ(a(t) -a(s), b), (voir, par exemple, [START_REF] Abdel-Hameed | A Gamma Wear Process[END_REF]). Voir [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF] Dans le but d'atténuer la dégradation du système au cours du temps et ainsi rallonger sa durée de vie, des actions de maintenance préventives peuvent être eectuées, en plus des maintenances correctives qui sont réalisées lors d'une défaillance. Du point de vue de la dégradation, un grand nombre de politiques de maintenance présentes dans la littérature considèrent des maintenances conditionnelles, où une maintenance préventive est eectuée dès lors que la dégradation atteint un seuil donné. Dans ce cas, la plupart des modèles de maintenance conditionnelle existants se limitent à des maintenances parfaites ("most of the existing CBM models have been limited to perfect maintenance actions"), comme relevé dans [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] (ou encore [START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF]). Cependant, plusieurs modèles de maintenance imparfaite apparaissent dans la littérature récente, dans le contexte de la dégradation de systèmes, voir [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF] pour une récente revue de la littérature. Certains modèles sont basés sur la notion d'âge virtuel, précédemment introduite dans le cadre des évènements récurrents (voir, par exemple, [START_REF] Giorgio | A new state-dependent degradation process and related model misidentication problems[END_REF][START_REF] Mercier | On the modelling of imperfect repairs for a continuously monitored gamma wear process through age reduction[END_REF]), pour lesquels le système est rajeuni par une maintenance. D'autres modèles supposent qu'un maintenance imparfaite réduit le niveau de dégradation du système, tels que [START_REF] Khatab | Condition-based selective maintenance for stochastically degrading multi-component systems under periodic inspection and imperfect maintenance[END_REF][START_REF] Letot | An adaptive degradation-based maintenance model taking into account both imperfect adjustments and agan replacements[END_REF][START_REF] Nicolai | Modelling and optimizing imperfect maintenance of coatings on steel structures[END_REF][START_REF] Ponchet | Maintenance policy on a nite time span for a gradually deteriorating system with imperfect improvements[END_REF][START_REF] Wu | A cost eective degradation-based maintenance strategy under imperfect repair[END_REF], et pouvant également être accompagnée par une augmentation de l'intensité de dégradation, comme dans [START_REF] Do | A proactive condition-based maintenance strategy with both perfect and imperfect maintenance actions[END_REF]. Dans d'autres articles, l'ecacité des maintenances est supposée décroissante avec le nombre de maintenances (voir, par exemple, [START_REF] Liu | An imperfect maintenance policy for mission-oriented systems subject to degradation and external shocks[END_REF][START_REF] Zhang | An ameliorated improvement factor model for imperfect maintenance and its goodness of t[END_REF]), et d'autres études plus poussées, comme dans [START_REF] Huynh | Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems[END_REF], considèrent des modèles de maintenance imparfaite tels que (i) les maintenances ont une ecacité aléatoire (ii) l'intensité de dégradation augmente avec le nombre de maintenance.

Dans tous ces articles, cependant, le point essentiel est l'optimisation de la politique de maintenance, en considérant à la fois des maintenances parfaites (remplacements) et imparfaites. À notre connaissance, peu d'articles traitent des modèles de maintenance imparfaite d'un point de vue statistique, excepté [START_REF] Zhang | Degradation-based maintenance decision using stochastic ltering for systems under imperfect maintenance[END_REF],

où les auteurs proposent un méthode du maximum de vraisemblance pour estimer les paramètres du processus de Wiener (décrivant la dégradation du système non maintenu) associée à une procédure itérative basée sur un ltre de Kalman pour les diérents facteurs impliqués dans les maintenances imparfaites successives.

Dans le cadre de systèmes se détériorant et sujets à des maintenances imparfaites, l'estimation des paramètres du processus de dégradation sous-jacent et de l'ecacité des maintenances est d'une grande utilité pour l'optimisation des politiques de maintenance. En eet, une fois les paramètres estimés, des prédictions sur le devenir du système maintenu peuvent être faites, permettant alors par exemple d'adapter (optimiser) la périodicité des actions de maintenance et de planier une révision générale.

D'un point de vue sécurité, la principale préoccupation est de s'assurer que les maintenances sont sufsamment ecaces, an de maintenir avec une forte probabilité le niveau de dégradation sous un seuil donné. Tant que ce seuil n'est pas atteint, les actions de maintenance peuvent être ajustées, soit en adaptant leur périodicité, soit en augmentant leur ecacité (lorsque c'est possible). Mis à part la sécurité, les coûts de maintenance est évidemment un autre axe important. C'est par exemple ce qui est fait dans [START_REF] Wu | A cost eective degradation-based maintenance strategy under imperfect repair[END_REF], où la minimisation des coûts est basée sur la durée de surveillance du système et de l'ecacité des actions de maintenance. Dans [START_REF] Huynh | Modeling past-dependent partial repairs for condition-based maintenance of continuously deteriorating systems[END_REF], l'auteur considère un seuil de dégradation au delà duquel un maintenance imparfaite est réalisée. L'optimisation est faite en fonction de ce seuil ainsi que la périodicité des inspections. Finalement, dans [START_REF] Wu | Maintenance policy for a system with a weighted linear combination of degradation processes[END_REF], une politique de maintenance est proposée, où un remplacement est eectué soit lorsque la dégradation dépasse un certain seuil, soit lorsqu'un nombre xé de maintenances imparfaites ont été réalisées. Voir, par exemple, ces trois articles et leurs références pour un aperçu en rapport avec l'optimisation de politiques de maintenance.

Cette thèse est dédiée au développement et à l'application de procédures d'estimation pour trois modèles spéciques de maintenance imparfaite dans le cadre de systèmes se dégradant suivant un processus gamma. Le document est divisé en quatre partie, incluant cette introduction.

Dans la Partie II, deux modèles de maintenance imparfaite sont étudiés. Il s'agit des modèles de réduction arithmétique de la dégradation d'ordre un et inni (ARD 1 et ARD ∞ , ARD signiant Arithmetic Reduction of Degradation), pour lesquels chaque maintenance réduit le niveau de dégradation du système. Le modèle ARD 1 a tout d'abord été introduit dans [START_REF] Castro | Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs[END_REF] et davantage étudié dans [34]. En se basant sur le modèle de réduction arithmétique de l'intensité (de défaillance) d'ordre un dans le cadre des évènements récurrents, l'idée de ce modèle est qu'une maintenance retire une proportion ρ de la dégradation accumulée par le système depuis la dernière maintenance (où ρ ∈ [0, 1[). Partant de cette même idée, [16] a également déni le modèle de réduction arithmétique de l'intensité d'ordre inni pour les évènements récurrents, que l'on adapte ici dans le contexte de la dégradation, conduisant à la dénition du modèle ARD ∞ . Concernant ce modèle, chaque maintenance réduit de ρ% le niveau de dégradation actuel, c'est-à-dire la dégradation accumulée par le système depuis sa mise en service. Une fois ces modèles dénis, on se place dans un cadre paramétrique et le schéma d'observation est lui aussi déni. La Il y a plusieurs points non abordés qu'il serait intéressant d'étudier en complément de cette thèse.

Il existe beaucoup d'autres modèles dans la continuité des modèles ARD étudiés ici, par exemple les modèles ARD m qui reposent sur l'idée qu'une action de maintenance réduit d'une certaine proportion la dégradation accumulée par le système depuis les m dernières maintenances. Une des axes de travail envisageable est de généraliser l'estimateur semi paramétrique de l'ecacité des maintenances à un tel modèle, ce qui semble possible car cette méthode repose sur la positivité des incréments du processus gamma. De plus, il serait intéressant d'adapter cette procédure d'estimation en considérant un processus de Lévy monotone autre que le processus gamma. Finalement, en se basant sur une rapide analyse numérique, il a été illustré que la condition qui assure une vitesse de convergence exponentielle est suf- Pour un modèle ARD 1 , une action de maintenance a pour eet de réduire d'une proportion donnée ρ ∈ [0, 1) la dégradation du système accumulée depuis la dernière maintenance. L'ecacité des actions de maintenance est alors mesurée par le paramètre ρ, ainsi, l'action de maintenance à l'instant jT réduit de ρ% la dégradation accumulée par le système durant l'intervalle de temps [ (j -1)T , jT [.

On note (Y t ) t≥0 le processus stochastique décrivant l'évolution du niveau de dégradation du système maintenu suivant le modèle ARD 1 . On note également (X (j) ) j∈N * une suite de copies indépendantes de (X t ) t≥0 , où X (j) correspond à la dégradation intrinsèque (sans maintenance) du système entre les instants (j -1)T et jT . Le système est supposé en parfait état à l'instant initial t = 0, soit Y 0 = X

(1) 0 = 0. Entre l'instant initial et la première maintenance à l'instant T , le système se dégrade suivant le processus X

T .

Pour t ∈ [T , 2T [, la dégradation du système correspond à la somme de la dégradation accumulée sur

T , et du niveau de dégradation juste après la première maintenance, donc

Finalement, pour tout t dans [nT , (n + 1)T [ et n dans N, on déduit une expression pour Y t , qui est : Ainsi, le niveau de dégradation Z t pour nT < t ≤ (n + 1)T , avec n dans N * , peut s'écrire comme

V ((n-1)T ) = Z (n-1)T + X En plus de cela, an de contourner le problème des intégrales en grande dimension, quatre méthodes alternatives basées sur le maximum de vraisemblance sont étudiées : la méthode du demi échantillon et le maximum de vraisemblance composite. La première méthode, basée sur les incréments, consiste à ne considérer qu'un incrément sur deux. Étant donné que le j-ème incrément est dépendant des j -1-ème et j + 1-ème, n'en choisir qu'un sur deux permet d'obtenir l'indépendance. Ainsi, suivant le principe du maximum de vraisemblance, une expression de la vraisemblance basée sur la moitié des incréments est donnée, nécessitant dans ce cas l'approximation d'intégrales de dimensions 1 seulement. Deux méthodes