# Importance and Sensitivity Analysis in Dynamic Reliability

Robert Eymard, Sophie Mercier, Michel Roussignol

► To cite this version:

Robert Eymard, Sophie Mercier, Michel Roussignol. Importance and Sensitivity Analysis in Dynamic Reliability. Methodology and Computing in Applied Probability, 2011, 13 (1), pp.75-104. 10.1007/s11009-009-9122-x. http://dx.doi.org/10.1007/s11009-009-9122-x

HAL Id: hal-00692992

https://hal.science/hal-00692992

Submitted on 31 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Importance and Sensitivity Analysis in Dynamic Reliability

Robert Eymard · Sophie Mercier · Michel Roussignol

**Abstract** In dynamic reliability, the evolution of a system is governed by a piece-wise deterministic Markov process, which is characterized by different input data. Assuming such data to depend on some parameter $p \in P$, our aim is to compute the first-order derivative with respect to each $p \in P$ of some functionals of the process, which may help to rank input data according to their relative importance, in view of sensitivity analysis. The functionals of interest are expected values of some function of the process, cumulated on some finite time interval $[0, t]$, and their asymptotic values per unit time. Typical quantities of interest hence are cumulated (production) availability, or mean number of failures on some finite time interval and similar asymptotic quantities. The computation of the first-order derivative with respect to $p \in P$ is made through a probabilistic counterpart of the adjoint state method, from the numerical analysis field. Examples are provided, showing the good efficiency of this method, especially in case of a large $P$.

**Keywords** Dynamic reliability · Sensitivity analysis · Importance factor · Piecewise deterministic Markov process · Chapman-Kolmogorov equation
# 1 Introduction

In dynamic reliability, the time-evolution of a system is described by a piecewise deterministic Markov process (PDMP) $(I_t, X_t)_{t \ge 0}$ introduced by Davis (1984). The first component $I_t$ is discrete, with values in a finite state space $E$. Typically, it indicates the state (up/down) for each component of the system at time $t$. The second component $X_t$, with values in a Borel set $V \subset \mathbb{R}^d$, stands for environmental conditions, such as temperature, pressure, and so on. Both components of the process interact one in each other: the process jumps at countably many isolated random times; by a jump from $(I_{t-}, X_{t-}) = (i, x)$ to $(I_t, X_t) = (j, y)$ (with $(i, x)$, $(j, y) \in E \times V$), the transition rate between the discrete states $i$ and $j$ depends on the environmental condition $x$ just before the jump and is a function $x \mapsto a(i, j, x)$. Similarly, the environmental condition just after the jump $X_t$ is distributed according to some distribution $\mu_{(i,j,x)}(dy)$, which depends on both components just before the jump $(i, x)$ and on the after jump discrete state $j$. Between jumps, the discrete component $I_t$ is constant, whereas the evolution of the environmental condition $X_t$ is deterministic, solution of a set of differential equations which depends on the fixed discrete state: given that $I_t = i$ for all $t \in [a, b]$, we have $\frac{d}{dt} X_t = v(i, X_t)$ for all $t \in [a, b]$, where $v$ is a mapping from $E \times V$ to $V$. Contrary to the general model from Davis (1993), we here assume that the eventual reaching of the frontier of $V$ does not entail jumps for the process $(I_t, X_t)_{t \ge 0}$.

Now, let $(I_t, X_t)_{t \ge 0}$ be a PDMP for which the jump rates $a(i, j, x)$, the jump distribution $\mu_{(i,j,x)}(dy)$ and the velocity fields $v(i, x)$ are assumed to depend on some family of parameters $P \in \mathbb{R}^k$, where $k \in \mathbb{N}$ can be quite large. Our aim here is to provide information about the sensitivity with respect to the elements of $P$, of expressions with the following form:

$$R_{\rho_0}(t) = \mathbb{E}_{\rho_0} \left( \int_0^t h(I_s, X_s) \, ds \right),$$

where $\rho_0$ is the initial distribution of the process and $h$ is some bounded measurable function which can also depend on $p \in P$. Such expressions include e.g. cumulative availability or production availability on some $[0, t]$, mean number of failures on $[0, t]$, mean time spent by $(X_s)_{0 \le s \le t}$ between two given bounds. We are also interested in the sensitivity with respect to $p \in P$ of the corresponding asymptotic quantities per unit time, namely in quantities of the form

$$\lim_{t \to +\infty} \frac{R_{\rho_0}(t)}{t}.$$

This sensitivity analysis can be guided by the knowledge of the first-order derivatives of $R_{\rho_0}(t)$ with respect to $p$, where $p \in P$. More specifically, we consider the following normalized derivative:

$$IF_p(t) = \frac{p}{R_{\rho_0}(t)} \frac{\partial R_{\rho_0}(t)}{\partial p}, \qquad (1)$$

for $t \le +\infty$ and $p \in P$, which we call the importance factor of parameter $p$ in $R_{\rho_0}(t)$. Such a dimensionless expression may help to rank the different parameters $p \in P$ according to their relative importance in $R_{\rho_0}(t)$. This kind of sensitivity analysis was already studied by Gandini (1990) and by Cao and Chen (1997) for pure jump
Markov processes with countable state space, and extended to PDMPs by Mercier and Roussignol (2008), with more restrictive a model than in the present paper however.

The starting point for this study is the Markov property of the process $(I_t, X_t)_{t \ge 0}$, which allows to write the associated Chapman–Kolmogorov equations fulfilled by its marginal distributions, see Davis (1993) or Cocozza-Thivent et al. (2006b). Such equations appear as weak forms of linear first order hyperbolic equations, see Eymard et al. (2008) e.g. The expressions for the derivatives of interest can then be obtained by solving the dual problem, as suggested by Lions (1968) for a wide class of partial differential equations. The Chapman–Kolmogorov equations and its dual problem are here solved using finite volume methods as provided by Cocozza-Thivent et al. (2006a) or Eymard et al. (2008), which here prove to be well adapted. An alternate way to compute the marginal distributions of the PDMP might be to start from the Markov renewal equations they fulfill, as proposed in Chiquet and Limnios (2008) for a specific case of PDMPs, or to use Monte-Carlo simulations.

In order to study the asymptotic quantities, one must put assumptions on the process to ensure its positive Harris recurrence. These assumptions may be proved to be true for Markov processes on general state space using techniques like those in Down et al. (1995) and Meyn and Tweedie (1993a, b), as is done for both examples at the end of this paper. Alternately, one may use specific results for piecewise deterministic Markov processes, as provided in Davis (1993) and Costa and Dufour (1999, 2003).

The paper is organized as follows: technical assumptions are set up in Section 2. An existence result for the derivatives of the transient quantities is provided in Section 3 and a computable expression is given for them in Section 4, using duality. Asymptotic quantities are studied in Section 5. Numerical procedures are exposed in Section 6, using finite volume methods. Two numerical studies close the paper in Sections 7 and 8.

## 2 Assumptions

All the paper is written under the following general assumptions:

- for all $i, j \in E$, the function $x \mapsto a(i, j, x)$ is continuous and bounded on $V$,

- for all $i, j \in E$ and all continuous and bounded function $f$ on $V$, the function $x \mapsto \int_V f(y) \mu_{(i,j,x)}(dy)$ is continuous and bounded on $V$,

- for all $i \in E$, the velocity field $v(i, x)$ is such that the function $x \mapsto v(i, x)$ is Lipschitz continuous and bounded by $V_1 = \max_{i \in E} \|v(i, \cdot)\|_\infty > 0$.

These assumptions guarantee the existence and uniqueness of the solution to the set of differential equations $\frac{dx}{dt} = v(i, x)$; we denote by $g(i, x, t)$ the single solution (defined on $E \times V \times \mathbb{R}$) such that $g(i, x, 0) = x$.

The jump rates $a(i, j, x)$, the jump distribution $\mu_{(i,j,x)}$, the velocity field $v(i, x)$ and the function $h$ are assumed to depend on some parameter $p$, where $p$ belongs to an open set $O \subset \mathbb{R}$ or $\mathbb{R}^k$. All the results are written in the case where $O \subset \mathbb{R}$ but extension to the case $O \subset \mathbb{R}^k$ is straightforward. We add exponent $(p)$ to each quantity depending on $p$, such as $h^{(p)}$ or $R_{\rho_0}^{(p)}(t)$.
Under the above technical assumptions, $(I_t^{(p)}, X_t^{(p)})_{t \ge 0}$ is known to be a Markov process with general state space $E \times V$, see Davis (1993), Cocozza-Thivent et al. (2006b), with strong Feller transition semi-group. We denote by $\rho_t^{(p)}(j, dy)$ the distribution of the process $(I_t^{(p)}, X_t^{(p)})_{t \ge 0}$ at time $t$ with initial distribution $\rho_0$ (independent of $p$) and by $P_t^{(p)}(i, x, j, dy)$ the transition probability distribution of $(I_t^{(p)}, X_t^{(p)})_{t \ge 0}$. We then have:

$$ 
\begin{align*}
R_{\rho_0}^{(p)}(t) &= \int_0^t \rho_s^{(p)} h^{(p)}(i, x) ds = \sum_{i \in E} \int_V \left( \int_0^t h^{(p)}(i, x) ds \right) \rho_s^{(p)}(i, dx) \\
&= \int_0^t \rho_0 P_s^{(p)} h^{(p)}(i, x) ds = \sum_{i \in E} \int_V \left( \int_0^t (P_s^{(p)} h^{(p)})(i, x) ds \right) \rho_0(i, dx)
\end{align*}
$$

In order to prove existence and to calculate derivatives of the functional $R_{\rho_0}^{(p)}$, we must give a sense to the derivatives of the transition probability distributions. With that aim, we need the following additional assumptions that we denote as **assumptions $H_1$ (resp. $H_2$):**

For each $p \in O$, there is some neighborhood $N(p)$ of $p$ in $O$ such that, for all $i, j \in E \times E$:

- the function $(x, p) \mapsto a^{(p)}(i, j, x)$ is bounded on $V \times N(p)$, once (resp. twice) continuously differentiable on $V \times O$, with all partial derivatives uniformly bounded on $V \times N(p)$,

- for all function $f^{(p)}(x)$ bounded and once (resp. twice) continuously differentiable on $V \times O$, with all partial derivatives uniformly bounded on $V \times N(p)$, the function $(x, p) \mapsto \int_V f^{(p)}(y) \mu_{(i,j,x)}^{(p)}(dy)$ is bounded and once (resp. twice) continuously differentiable on $V \times O$, with all partial derivatives uniformly bounded on $V \times N(p)$,

- the function $(x, p) \mapsto v^{(p)}(i, x)$ is bounded on $V \times N(p)$, once (resp. twice) continuously differentiable on $V \times O$, with all partial derivatives uniformly bounded on $V \times N(p)$,

- the function $(x, p) \mapsto h^{(p)}(i, x)$ is bounded on $V \times N(p)$, once (resp. twice) continuously differentiable on $V \times O$ with uniformly bounded partial derivatives on $V \times N(p)$.

The third point implies (see e.g. Cartan 1967) that, for all $i \in E$, the function $(x, p) \mapsto g^{(p)}(i, x, s)$, solution of the differential equation $\frac{dg}{dt} = v^{(p)}(i, g)$ with $g(i, x, 0) = x$, is once (resp. twice) continuously differentiable on $V \times N(p)$, with all partial derivatives uniformly bounded on $V \times N(p)$.

Throughout the paper, under assumptions $H_1$ or $H_2$, for each $p \in O$, we shall refer to a $N(p)$ fulfilling the four points of the assumptions without any further notice.

We may now give a sense to the derivatives of the transition probability distributions, which is done in next section.
### 3 Existence Result

We shall use the infinitesimal generators of both Markov processes $(I_t, X_t)_{t \ge 0}$ and $(I_t, X_t, t)_{t \ge 0}$:

**Definition 1** Let $\mathcal{D}_{H_0}$ be the set of functions $f(i, x)$ from $E \times V$ to $\mathbb{R}$ such that for all $i \in E$ the function $x \mapsto f(i, x)$ is bounded, continuously differentiable on $V$ and such that the function $x \mapsto \mathbf{v}^{(p)}(i, x) \cdot \nabla f(i, x)$ is bounded on $V$.

For $f \in \mathcal{D}_{H_0}$, we define

$$H_0^{(p)} f(i, x) = \sum_{j \in E} a^{(p)}(i, j, x) \int_V f(j, y) \mu_{(i,j,x)}^{(p)}(dy) + \mathbf{v}^{(p)}(i, x) \cdot \nabla f(i, x)$$

where we set $a^{(p)}(i, i, x) = - \sum_{j \neq i} a^{(p)}(i, j, x)$ and $\mu_{(i,i,x)}^{(p)} = \delta_x$.

Let $\mathcal{D}_H$ be the set of functions $f(i, x, s)$ from $E \times V \times \mathbb{R}_+$ to $\mathbb{R}$ such that for all $i \in E$ the function $(x, s) \mapsto f(i, x, s)$ is bounded, continuously differentiable on $V \times \mathbb{R}_+$ and such that the function $(x, s) \mapsto \frac{\partial f}{\partial s}(i, x, s) + \mathbf{v}^{(p)}(i, x) \cdot \nabla f(i, x, s)$ is bounded on $V \times \mathbb{R}_+$.

For $f \in \mathcal{D}_H$, we define

$$H^{(p)} f(i, x, s) = \sum_{j} a^{(p)}(i, j, x) \int_{V} f(j, y, s) \mu_{(i,j,x)}^{(p)}(dy) + \frac{\partial f}{\partial s}(i, x, s) \\ + \mathbf{v}^{(p)}(i, x) \cdot \nabla f(i, x, s) \quad (2)$$

We then have for all $f \in \mathcal{D}_{H_0}$:

$$P_t^{(p)} f = f + \int_0^t H_0^{(p)} P_u^{(p)} f \, du \quad (3)$$

$$\rho_t^{(p)} f = \rho_0 f + \int_0^t \rho_u^{(p)} H_0^{(p)} f \, du \quad (4)$$

and for all $f \in \mathcal{D}_H$:

$$P_t^{(p)} f(\cdot, \cdot, t) = f(\cdot, \cdot, 0) + \int_0^t P_u^{(p)} H^{(p)} f(\cdot, \cdot, u) \, du \quad (5)$$

$$\rho_t^{(p)} f(\cdot, \cdot, t) = \rho_0 f(\cdot, \cdot, 0) + \int_0^t \rho_u^{(p)} H^{(p)} f(\cdot, \cdot, u) \, du \quad (6)$$

These are Chapman-Kolmogorov equations.

Thanks to these equations and to Theorem 4 in the Appendix, we get the following result:

**Proposition 1** Let $f(i, x)$ be a function (independent of $p$) from $E \times V$ to $\mathbb{R}$ such that for all $i \in E$, the function $x \mapsto f(i, x)$ is bounded, continuously differentiable on $V$ with uniformly bounded partial derivatives on $V$. Under assumptions $\mathcal{H}_1$, for all $i \in E$, the function $(x, p) \mapsto (P_t^{(p)} f)(i, x)$ is continuously differentiable on $V \times O$, with all partial derivatives uniformly bounded on $V \times N(p)$ for all $p \in O$. Under assumptions
$\mathcal{H}_2$, if for all $i \in E$ the function $x \mapsto f(i, x)$ is bounded, twice continuously differentiable on $V$ with uniformly bounded partial derivatives on $V$, then for all $i \in E$, the function $(x, p) \mapsto (P_t^{(p)} f)(i, x)$ is twice continuously differentiable on $V \times O$, with all partial derivatives uniformly bounded on $V \times N(p)$ for all $p \in O$.

**Proof** Let $p_0 \in P$ and $N(p_0)$ be a neighborhood of $p_0$. Setting $\varphi(s, i, (x, p)) = P_t^{(p)} f(i, x)$ for $(i, x, p, s) \in E \times V \times N(p_0) \times \mathbb{R}_+$, the Chapman-Kolmogorov Eq. 3 can be written as:

$$ \varphi(t, i, (x, p)) = f(i, x) + \int_0^t \left( \sum_{j \in E} a^{(p)}(i, j, x) \int_V \varphi(u, j, (y, p)) \mu_{(i,j,x)}^{(p)}(dy) + \mathbf{v}^{(p)}(i, x) \cdot \nabla_x \varphi(u, i, (x, p)) \right) du \quad (7) $$

where $\nabla_x$ stands for the gradient with respect of $x$. Using a similar method as in Cocozza-Thivent et al. (2006b), we introduce the function $\tilde{\varphi}$ defined by:

$$ \tilde{\varphi}(t, i, (x, p)) = \varphi(t, i, (g^{(p)}(i, x, -t), p)) $$

Noting that

$$ \varphi(t, i, (x, p)) = \tilde{\varphi}(t, i, (g^{(p)}(i, x, t), p)), $$

Equation 7 may now be written as:

$$ \begin{aligned} & \tilde{\varphi}(t, i, (x, p)) \\ & = f(i, x) \\ & + \int_0^t \left( \sum_{j \in E} a^{(p)}(i, j, g^{(p)}(i, x, -u)) \int_V \tilde{\varphi}(u, j, (g^{(p)}(j, y, u), p)) \mu_{(i,j,g^{(p)}(i,x,-u))}^{(p)}(dy) \right) du \end{aligned} $$

Using the notations of Theorem 4, we set $z = (x, p)$, $\tilde{\varphi}_0(i, z) = f(i, x)$ and

$$ \mathcal{F}(s, \psi)(i, z) = \sum_{j \in E} a^{(p)}(i, j, g^{(p)}(i, x, -u)) \int_V \psi(j, (g^{(p)}(j, y, s), p)) \mu_{(i,j,g^{(p)}(i,x,-s))}^{(p)}(dy) $$

for $(i, x, p, s) \in E \times V \times N(p_0) \times \mathbb{R}_+$ and $\psi: E \times V \times N(p_0) \to \mathbb{R}$ bounded and measurable. The function $\tilde{\varphi}$ then satisfies the equations:

$$ \tilde{\varphi}(t, i, z) = \tilde{\varphi}_0(i, z) + \int_0^t \mathcal{F}(u, \tilde{\varphi}(u, \cdot, \cdot))(i, z) du. $$

We now check that assumptions of Theorem 4 are fulfilled and we set $I$ to be an interval such that $[0, T] \subset I$. Symbols $W = V \times N(p_0)$, $C_1^b(E, W)$, $C_1^b(I, E, W)$
and $J(\mathcal{F})$ are as in the Appendix. With such notations, the function $\tilde{\phi}_0$ belongs to $C_1^b(E, W)$ and assumptions $\mathcal{H}_1$ imply that:

- if $\Psi \in C_1^b(I, E, W)$, then $\mathcal{F}(s, \Psi(s, \cdot, \cdot))(i, (x, p))$ is in $C_1^b(I, E, W)$,

- if $\Psi \in C_1^b(E, W)$, then for all $s \in I$, $\mathcal{F}(s, \Psi)$ and $\nabla \mathcal{F}(s, \Psi)$ are uniformly Lipschitz with respect to $\Psi$.

Assumptions $\mathcal{H}_2$ and $\tilde{\phi}_0 \in C_2^b(E, W)$ imply that:

- if $\Psi \in C_2^b(I, E, W)$, then $\mathcal{F}(s, \Psi(s, \cdot, \cdot))(i, (x, p))$ is in $C_2^b(I, E, W)$,

- if $\Psi \in C_2^b(E, W)$ then for all $s \in I$, $\mathcal{F}(s, \Psi)$, $\nabla \mathcal{F}(s, \Psi)$ and $J(\mathcal{F})(s, \Psi)$ are uniformly Lipschitz with respect to $\Psi$.

All required assumptions are then checked, which provides the result. $\square$

**Remark 1** Using the explicit form of transition probabilities (see Cocozza-Thivent et al. 2006b), the functional $f \mapsto \frac{\partial}{\partial p}(P_t^{(p)} f)(i, x)$ appears to be a continuous linear functional of $f$ which involves $f$ and $\frac{\partial f}{\partial x_k}$ for $1 \le k \le d$ (and possibly other derivatives). Limiting this functional to the set of functions $f$ such that $x \mapsto f(i, x)$ is infinitely differentiable with a compact support for all $i \in E$, it may then be seen as a distribution. We actually define this functional on the greater space of bounded continuously differentiable functions with uniformly bounded partial derivatives and we use the following notation:

$$ \frac{\partial}{\partial p} (P_t^{(p)} f)(i, x) = \left\langle \frac{\partial P_t^{(p)}}{\partial p}, f \right\rangle (i, x) $$

Next corollary is straightforward.

**Corollary 1** Under assumptions $\mathcal{H}_1$, the function $(x, p) \mapsto (P_t^{(p)} h^{(p)})(i, x)$ is continuously differentiable on $V \times O$ for all $i \in E$, with all partial derivatives uniformly bounded on $V \times N(p)$ for all $p \in O$. Besides:

$$ \frac{\partial}{\partial p} (P_t^{(p)} h^{(p)})(i, x) = \left\langle \frac{\partial P_t^{(p)}}{\partial p}, h^{(p)} \right\rangle (i, x) + \left( P_t^{(p)} \frac{\partial h^{(p)}}{\partial p} \right) (i, x) $$

for all $(i, x, t, p) \in E \times V \times \mathbb{R}^+ \times O$.
We derive the following theorem:

**Theorem 1** Under assumptions $\mathcal{H}_1$, the function $p \mapsto R_{\rho_0}^{(p)}(t)$ is continuously differentiable with respect to $p$ and

$$ 
\begin{aligned}
\frac{\partial}{\partial p} (R_{\rho_0}^{(p)}(t)) &= \int_0^t \sum_{i \in E} \int_V \left\langle \frac{\partial P_u^{(p)}}{\partial p}, h^{(p)} \right\rangle (i, x) \rho_0(i, dx) \, du \\
&\quad + \int_0^t \sum_{i \in E} \sum_{j \in E} \int_V \int_V \frac{\partial}{\partial p} h^{(p)}(j, y) P_u^{(p)}(i, x, j, dy) \rho_0(i, dx) \, du \\
&= \int_0^t \left\langle \frac{\partial \rho_u^{(p)}}{\partial p}, h^{(p)} \right\rangle \, du + \int_0^t \rho_u^{(p)} \frac{\partial h^{(p)}}{\partial p} \, du
\end{aligned}
\quad (8) $$

where we set:

$$ \left\langle \frac{\partial \rho_u^{(p)}}{\partial p}, h^{(p)} \right\rangle = \sum_{i \in E} \int_V \left\langle \frac{\partial P_u^{(p)}}{\partial p}, h^{(p)} \right\rangle (i, x) \rho_0(i, dx). \quad (9) $$

*Remark 2* Assumption on $h^{(p)}$ in $\mathcal{H}_1$ may be generalized to the case where $h^{(p)}(i, x)$ is continuously differentiable only for almost all (a.a.) $x \in V$, where a.a. means with respect to Lebesgue measure.

Our purpose now is to compute this derivative. The marginal distribution $\rho_u^{(p)}(j, dy)$ may be estimated by different methods, as indicated in the introduction. However, we do not know how to compute directly the derivative of the marginal distribution which appear in the above expression. We now transform it in order to make it easier to compute.

## 4 Results in the Transient Case

Let us first define the notion of importance functions.

**Definition 2** We say that a function $\varphi_t^{(p)} \in \mathcal{D}_H$ is the importance function associated to the function $h^{(p)}$ and to $t$ if:

- $\varphi_t^{(p)}$ is solution of the differential equation $H^{(p)}\varphi_t^{(p)}(i, x, s) = h^{(p)}(i, x)$ for all $(i, x, s) \in E \times V \times [0, t[$,

- $\varphi_t^{(p)}(i, x, t) = 0$ for all $(i, x) \in E \times V$.

**Proposition 2** Let $t > 0$. Under assumptions $\mathcal{H}_1$, the importance function associated to $(h^{(p)}, t)$ exists, is unique, is continuously differentiable on $V \times O$, bounded with all partial derivatives uniformly bounded on $V \times N(p)$ for all $p \in O$. Under assumptions $\mathcal{H}_2$, the importance function associated to $(h^{(p)}, t)$ is moreover twice continuously differentiable on $V \times O$ with all partial derivatives uniformly bounded on $V \times N(p)$ for all $p \in O$.
Proof Let $\varphi^{(p)} \in \mathcal{D}_H$ (eventually depending on $p$) and let us set:

$$ \tilde{\varphi}^{(p)}(i, z, s) = \varphi^{(p)}(i, g(i, z, -s), t-s) \quad (10) $$

for all $(i, z, s) \in E \times V \times [0, t]$. The function $\tilde{\varphi}^{(p)}$ is bounded, continuously differentiable on $V$ with partial derivative with respect of time given by:

$$ \frac{\partial \tilde{\varphi}^{(p)}}{\partial s_1}(i, z, s) = -\frac{\partial \varphi^{(p)}}{\partial s_1}(i, g^{(p)}(i, z, -s), t-s) \\ -\nabla \varphi^{(p)}(i, g^{(p)}(i, z, -s), t-s) \cdot \mathbf{v}^{(p)}(i, g^{(p)}(i, z, -s)) \quad (11) $$

The function $\varphi^{(p)}$ and $\tilde{\varphi}^{(p)}$ are in one-to-one correspondence with:

$$ \varphi^{(p)}(i, z, t-s) = \tilde{\varphi}^{(p)}(i, g^{(p)}(i, z, s), s) \quad (12) $$

for all $(i, z, s) \in E \times V \times [0, t]$. Using Eqs. 10–11, we get:

$$ \begin{aligned} & H^{(p)} \varphi^{(p)}(i, g^{(p)}(i, z, -s), t-s) \\ &= \sum_j a^{(p)}(i, j, g^{(p)}(i, z, -s)) \int_V \varphi^{(p)}(j, y, t-s) \mu^{(p)}_{(i,j,g^{(p)}(i,z,-s))}(dy) - \frac{\partial \tilde{\varphi}^{(p)}}{\partial s_1}(i, z, s) \\ &= \tilde{H}^{(p)} \tilde{\varphi}^{(p)}(i, z, s) \end{aligned} $$

where

$$ \begin{aligned} & \tilde{H}^{(p)} \tilde{\varphi}^{(p)}(i, z, s) \\ &= \sum_{j \in E} a^{(p)}(i, j, g^{(p)}(i, z, -s)) \int_V \tilde{\varphi}^{(p)}(j, g^{(p)}(j, y, s), s) \mu^{(p)}_{(i,j,g^{(p)}(i,z,-s))}(dy) \\ & \quad - \frac{\partial \tilde{\varphi}^{(p)}}{\partial s_1}(i, z, s). \end{aligned} $$

The problem now resumes to show existence, uniqueness and regularity of $\tilde{\varphi}^{(p)}$ such that

$$ \tilde{H}^{(p)} \tilde{\varphi}^{(p)}(i, z, s) = h^{(p)}(i, g^{(p)}(i, z, -s)) \quad (13) $$

for all $(i, z, s) \in E \times V \times [0, t]$ with $\tilde{\varphi}^{(p)}(i, z, 0) = \varphi^{(p)}(i, z, t) = 0$ for all $(i, z) \in E \times V$. The o.d.e. Eq. 13 may be written as:

$$ \tilde{\varphi}^{(p)}(i, z, t) = \int_{0}^{t} \mathcal{R}^{(p)}(\tilde{\varphi}^{(p)}(\cdot, \cdot, s), s)(i, z) \, ds $$

with

$$ \begin{aligned} & \mathcal{R}^{(p)}(\tilde{\varphi}^{(p)}(\cdot, \cdot, s), s)(i, z) \\ &= \sum_{j \in E} a^{(p)}(i, j, g^{(p)}(i, z, -s)) \int_V \tilde{\varphi}^{(p)}(j, g^{(p)}(j, y, s), s) \mu^{(p)}_{(i,j,g^{(p)}(i,z,-s))}(dy) \\ & \quad - h^{(p)}(i, g^{(p)}(i, z, -s)). \end{aligned} $$

Theorem 4 from the Appendix then provides the result in the same way as for Proposition 1. □
Next lemma provides a duality result based on the definition of $\varphi_t^{(p)}$, which transports the differentiation with respect of $p$ from the marginal distribution $\rho_s^{(p)}$ to operator $H^{(p)}$.

**Lemma 1** Under assumptions $\mathcal{H}_2$, we have:

$$ \int_0^t \left\langle \frac{\partial \rho_s^{(p)}}{\partial p}, H^{(p)} \varphi_t^{(p)}(\cdot, \cdot, s) \right\rangle ds = - \int_0^t \rho_s^{(p)} \frac{\partial H^{(p)}}{\partial p} \varphi_t^{(p)}(\cdot, \cdot, s) ds $$

where we set:

$$ \begin{aligned} \frac{\partial H^{(p)}}{\partial p} \varphi (i, x, s) := & \sum_{j \in E} \frac{\partial a^{(p)}}{\partial p}(i, j, x) \int_V \varphi(j, y, s) \mu_{(i,j,x)}^{(p)}(dy) \\ & + \sum_{j \in E} a^{(p)}(i, j, x) \frac{\partial}{\partial p} \left( \int_V \varphi(j, y, s) \mu_{(i,j,x)}^{(p)}(dy) \right) \\ & + \frac{\partial v^{(p)}}{\partial p}(i, x) \cdot \nabla \varphi(i, x, s) \end{aligned} $$

for all $\varphi \in \mathcal{D}_H$ and all $(i, x, s) \in E \times V \times \mathbb{R}_+$.

*Proof* Thanks to Proposition 2, functions $\varphi_t^{(p)}(i, x, s)$ and $\frac{\partial \varphi_t^{(p)}}{\partial p}(i, x, s)$ are in $\mathcal{D}_H$. We first know from the Chapman-Kolmogorov Eq. 6 applied to $\varphi_t^{(p)}$ that:

$$ \int_0^t \rho_s^{(p)} H^{(p)} \varphi_t^{(p)}(\cdot, \cdot, s) ds = \rho_t^{(p)} \varphi_t^{(p)}(\cdot, \cdot, t) - \rho_0 \varphi_t^{(p)}(\cdot, \cdot, 0) = -\rho_0 \varphi_t^{(p)}(\cdot, \cdot, 0) $$

due to $\varphi_t^{(p)}(\cdot, \cdot, t) = 0$. By differentiating this expression with respect to $p$, we derive:

$$ \int_0^t \left\langle \frac{\partial \rho_s^{(p)}}{\partial p}, H^{(p)} \varphi_t^{(p)}(\cdot, \cdot, s) \right\rangle ds + \int_0^t \rho_s^{(p)} \frac{\partial H^{(p)}}{\partial p} \varphi_t^{(p)}(\cdot, \cdot, s) ds \\ + \int_0^t \rho_s^{(p)} H^{(p)} \frac{\partial \varphi_t^{(p)}}{\partial p}(\cdot, \cdot, s) ds = -\rho_0 \frac{\partial \varphi_t^{(p)}}{\partial p}(\cdot, \cdot, 0). \quad (14) $$

Chapman-Kolmogorov Eq. 6 applied to $\frac{\partial \varphi_t^{(p)}}{\partial p}$ gives:

$$ \begin{aligned} \int_0^t \rho_s^{(p)} H^{(p)} \frac{\partial \varphi_t^{(p)}}{\partial p}(\cdot, \cdot, s) ds &= \rho_t^{(p)} \frac{\partial \varphi_t^{(p)}}{\partial p}(\cdot, \cdot, t) - \rho_0 \frac{\partial \varphi_t^{(p)}}{\partial p}(\cdot, \cdot, 0) \\ &= -\rho_0 \frac{\partial \varphi_t^{(p)}}{\partial p}(\cdot, \cdot, 0). \end{aligned} $$

Hence the result, substituting this last outcome in Eq. 14. $\square$

We easily derive the following result, using the definition of $R_{\rho_0}^{(p)}(t)$ and of $\varphi_t^{(p)}(\cdot, \cdot, s)$.
**Theorem 2** Under assumptions $\mathcal{H}_2$, we have:

$$ \frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t) = \int_0^t \rho_s^{(p)} \frac{\partial h^{(p)}}{\partial p} ds - \int_0^t \rho_s^{(p)} \frac{\partial H^{(p)}}{\partial p} \varphi_t^{(p)} (\cdot, \cdot, s) ds, \quad (15) $$

where $\varphi_t^{(p)}$ is the importance function associated to $(h^{(p)}, t)$.

**Remark 3** In case that $R_{\rho_0}(t)$ depends on a family of parameters $P$, the computation of all $\frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t)$ for each $p \in P$ consequently requires only one single computation of both $\rho_s^{(p)}$ and $\varphi_t^{(p)}$, and simple summations on $[0, t]$ for each value of $p$. This is to be compared to the usual finite differences method, which requires one first evaluation of $\rho_s^{(p)}$ and another evaluation of $\rho_s^{(p, \epsilon)}$ for each value of $p \in P$ and some $\epsilon > 0$, where $P_{p, \epsilon}$ stands for the family $P$ in which parameter $p$ has been changed into $p + \epsilon$. In case of a large $P$, the present methods, which only requires two effective computations, consequently appears as cheaper than finite differences.

In applications, the importance function will generally be computed numerically. An analytical form is however available, which is also useful for the asymptotic study:

**Lemma 2** Let $t > 0$. Under assumptions $\mathcal{H}_1$, the importance function associated to $(h^{(p)}, t)$ is given by:

$$ \varphi_t^{(p)}(i, x, s) = \begin{cases} -\int_0^{t-s} (P_u^{(p)} h^{(p)})(i, x) du & \text{if } 0 \le s \le t \\ 0 & \text{otherwise} \end{cases} \quad (16) $$

for all $(i, x) \in E \times V$.

*Proof* Let $\varphi_t^{(p)}$ be the function defined by Eq. 16. It is clear that $\varphi_t^{(p)}(i, x, t) = 0$ for all $(i, x) \in E \times V$ and thanks to Proposition 1, the function $\varphi_t^{(p)}$ is in $\mathcal{D}_H$. Besides, for $0 \le s < t$, we have:

$$ \begin{aligned} (H^{(p)} \varphi_t^{(p)})(\cdot, \cdot, s) &= -H^{(p)} \int_0^{t-s} P_u^{(p)} h^{(p)} du \\ &= -H_0^{(p)} \left( \int_0^{t-s} P_u^{(p)} h^{(p)} du \right) - \frac{\partial}{\partial s} \left( \int_0^{t-s} P_u^{(p)} h^{(p)} du \right) \\ &= -\int_0^{t-s} H_0^{(p)} (P_u^{(p)} h^{(p)}) du + P_{t-s}^{(p)} h^{(p)} \\ &= h^{(p)} \end{aligned} $$

due to the Chapman-Kolmogorov Eq. 3, which completes the proof. $\square$

We easily derive the following Corollary from Theorem 2 and the previous lemma.

**Corollary 2** Under assumptions $\mathcal{H}_2$, we have:

$$ \frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t) = \int_0^t \rho_s^{(p)} \frac{\partial h^{(p)}}{\partial p} ds + \int_0^t \rho_s^{(p)} \frac{\partial H^{(p)}}{\partial p} \left( \int_0^{t-s} P_u^{(p)} h^{(p)} du \right) ds \quad (17) $$
Equation 17 is an extension of the results of Gandini (1990) for pure jump Markov processes with countable state space.

## 5 Asymptotic Results

We are now interested in asymptotic results and we need to assume the process $(I_t, X_t)_{t \ge 0}$ to be uniformly ergodic, according to the following assumptions $\mathcal{H}_3$:

- the process $(I_t, X_t)_{t \ge 0}$ is positive Harris-recurrent with $\pi^{(p)}$ as unique stationary distribution,

- for each $p \in O$, there exists a function $f^{(p)}$ such that

$$ \int_0^{+\infty} f^{(p)}(u) \, du < +\infty, \quad \int_0^{+\infty} u f^{(p)}(u) \, du < +\infty, \quad \lim_{u \to +\infty} f^{(p)}(u) = 0, \quad (18) $$

and

$$ \left| (P_u^{(p)} h^{(p)}) (i, x) - \pi^{(p)} h^{(p)} \right| \le f^{(p)}(u) \quad (19) $$

for all $(i, x) \in E \times V$ and all $u \ge 0$.

As already noted in the introduction, such assumptions may be proved to be true by using results for Markov processes on general state spaces or specific results for piecewise deterministic Markov processes (see the introduction for references).

In order not to give too technical assumptions difficult to check in practice, we constrain our asymptotic study to the special case where only the jump rates $a^{(p)}(i, j, x)$ and the function $h^{(p)}(i, x)$ depend on the parameter $p$. The quantities $\mu_{(i,j,x)}$ and $v(i, x)$ are consequently assumed to be independent of $p$. Assumptions $\mathcal{H}_2$ are then substituted by assumptions $\mathcal{H}_2'$, where conditions on $\mu_{(i,j,x)}$ and on $v(i, x)$ are removed.

We now transform Eq. 17 in view of studying its asymptotic expression.

**Lemma 3** Under assumptions $\mathcal{H}_2$ and $\mathcal{H}_3$, we have:

$$ \frac{1}{t} \frac{\partial R_{\rho_0}^{(p)}}{\partial p} (t) = \frac{1}{t} \int_0^t \rho_s^{(p)} \frac{\partial h^{(p)}}{\partial p} ds \\ + \frac{1}{t} \int_0^t \rho_s^{(p)} \frac{\partial H^{(p)}}{\partial p} \left( \int_0^{t-s} \left( P_u^{(p)} h^{(p)} - \pi^{(p)} h^{(p)} \right) du \right) ds \quad (20) $$

*Proof* The first term is clear. Besides, setting **1** to be the constant function equal to 1, we have: $\frac{\partial H^{(p)}}{\partial p} \mathbf{1} = 0$ since $H^{(p)} \mathbf{1} = 0$. As $\pi^{(p)} h^{(p)}$ is a constant (independent of $(i, x)$), we derive

$$ \frac{\partial H^{(p)}}{\partial p} (\pi^{(p)} h^{(p)}) = (\pi^{(p)} h^{(p)}) \frac{\partial H^{(p)}}{\partial p} \mathbf{1} = 0 $$

and consequently:

$$ \frac{\partial H^{(p)}}{\partial p} \left( \int_{0}^{t-s} \pi^{(p)} h^{(p)} du \right) = (t-s) \frac{\partial H^{(p)}}{\partial p} (\pi^{(p)} h^{(p)}) = 0. $$

Whence the result. □
We may now prove existence and provide an asymptotic expression for $\frac{1}{t} \frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t)$.

**Theorem 3** Let us assume that $\mu_{(i,j,x)}$ and $v(i, x)$ are independent of $p$ and that $H'_2$, $H_3$ are true. Then:

$$ U h^{(p)}(i, x) := \int_{0}^{+\infty} \left( (P_{u}^{(p)} h^{(p)})(i, x) - \pi^{(p)} h^{(p)} \right) du $$

exists for all $(i, x) \in E \times V$ and $U h^{(p)}$ is element of $D_{H_0}$. Besides:

$$ \lim_{t \to +\infty} \frac{1}{t} \frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t) = \pi^{(p)} \frac{\partial h^{(p)}}{\partial p} + \pi^{(p)} \frac{\partial H_0^{(p)}}{\partial p} U h^{(p)} \quad (21) $$

where we set:

$$ \frac{\partial H_0^{(p)}}{\partial p} \varphi_0(i, x) := \sum_{j \in E} \frac{\partial a^{(p)}}{\partial p}(i, j, x) \int_V \varphi_0(j, y) \mu_{(i,j,x)}(dy) $$

for all $\varphi_0 \in D_{H_0}$ and all $(i, x) \in E \times V$.

*Proof* Due to $H_3$, the quantity $Uh^{(p)}(i, x)$ is clearly defined for all $(i, x) \in E \times V$ and $Uh^{(p)}$ is element of $D_{H_0}$. Under $H_3$, we know that:

$$ \lim_{t \to +\infty} \frac{1}{t} \int_{0}^{t} \rho_{s}^{(p)} \varphi^{(p)} ds = \pi^{(p)} \varphi^{(p)} \quad (22) $$

for each measurable and bounded $\varphi^{(p)}$. The first term in right side of Eq. 20 consequently converges to the first term in Eq. 21. For the second term, setting $U_s := \int_s^{+\infty} (P_u^{(p)} h^{(p)} - \pi^{(p)} h^{(p)}) du$, we have:

$$ \frac{1}{t} \int_{0}^{t} \rho_{s}^{(p)} \frac{\partial H^{(p)}}{\partial p} \left( \int_{0}^{t-s} (P_{u}^{(p)} h^{(p)} - \pi^{(p)} h^{(p)}) du \right) ds \\ = \frac{1}{t} \int_{0}^{t} \rho_{s}^{(p)} \frac{\partial H^{(p)}}{\partial p} U h^{(p)} ds - \frac{1}{t} \int_{0}^{t} \rho_{s}^{(p)} \frac{\partial H^{(p)}}{\partial p} U_{t-s} ds. $$

By assumption, the function $\frac{\partial H^{(p)}}{\partial p} Uh^{(p)}$ is bounded and independent of time, and consequently:

$$ \lim_{t \to +\infty} \frac{1}{t} \int_{0}^{t} \rho_{s}^{(p)} \frac{\partial H^{(p)}}{\partial p} U h^{(p)} ds = \pi^{(p)} \frac{\partial H^{(p)}}{\partial p} U h^{(p)} = \pi^{(p)} \frac{\partial H_{0}^{(p)}}{\partial p} U h^{(p)}. $$

It now remains to prove that

$$ \lim_{t \to +\infty} \frac{1}{t} \int_{0}^{t} \rho_{s}^{(p)} \frac{\partial H^{(p)}}{\partial p} U_{t-s} ds = 0. $$

We have:

$$ \left| \frac{\partial H^{(p)}}{\partial p} U_{t-s}(i, x) \right| \le \sum_{j \in E} \left| \frac{\partial a^{(p)}}{\partial p}(i, j, x) \right| \left| \int_V U_{t-s}(j, y) \mu_{(i,j,x)}(dy) \right| $$
As $|U_s(i, x)| \le \int_s^{+\infty} f^{(p)}(u) du$ due to $\mathcal{H}_3$, we have:

$$ \left| \frac{\partial H^{(p)}}{\partial p} U_{t-s}(i, x) \right| \le K \int_{t-s}^{+\infty} f^{(p)}(u) du $$

and

$$ \begin{aligned} \left| \frac{1}{t} \int_0^t \rho_s^{(p)} \frac{\partial H^{(p)}}{\partial p} U_{t-s} ds \right| &\le K \frac{1}{t} \int_0^t \left( \int_{t-s}^{+\infty} f^{(p)}(u) du \right) ds \\ &\le K \frac{1}{t} \int_0^{+\infty} u f^{(p)}(u) du. \end{aligned} $$

This ends the proof.

□

The previous theorem provides an extension of the results given in Cao and Chen (1997) for pure jump Markov processes with countable state space.

Next proposition gives a tool to compute the function $Uh^{(p)}$.

**Proposition 3** Let us assume that $\mu_{(i,j,x)}$ and $v(i, x)$ are independent of $p$ and that $\mathcal{H}_2$, $\mathcal{H}_3$ are true. Then the function $Uh^{(p)}$ is the unique solution of the following ordinary differential equation:

$$ H_0^{(p)} Uh^{(p)}(i, x) = \pi^{(p)} h^{(p)} - h^{(p)}(i, x) \quad (23) $$

for all $(i, x) \in E \times V$ such that $\pi^{(p)}Uh^{(p)} = 0$.

*Proof* Under assumptions $\mathcal{H}_3$, we clearly have $\pi^{(p)}Uh^{(p)} = 0$ by Fubini's theorem. We now check that $Uh^{(p)}$ is solution of Eq. 23: under $\mathcal{H}_3$, we may write

$$ \begin{aligned} H_0^{(p)} Uh^{(p)} (i, x) &= H_0^{(p)} \int_0^{+\infty} \left( (P_u^{(p)} h^{(p)})(i, x) - \pi^{(p)} h^{(p)} \right) du \\ &= \int_0^{+\infty} H_0^{(p)} P_u^{(p)} h^{(p)} (i, x) du \end{aligned} $$

since $H_0^{(p)}(\pi^{(p)}h^{(p)}) = (\pi^{(p)}h^{(p)}) H_0^{(p)} \mathbf{1} = 0$. We derive:

$$ \begin{aligned} H_0^{(p)} Uh^{(p)} (i, x) &= \lim_{t \to +\infty} \int_0^t (H_0^{(p)} P_u^{(p)} h^{(p)} (i, x)) du \\ &= \lim_{t \to +\infty} P_t^{(p)} h^{(p)} (i, x) - h^{(p)} (i, x) \text{ (Chapman-Kolmogorov equation)} \\ &= \pi^{(p)} h^{(p)} - h^{(p)} (i, x), \end{aligned} $$

which shows that $Uh^{(p)}$ is solution of Eq. 23. Now, let $\varphi$ be the difference between two solutions $\varphi_1, \varphi_2$ of Eq. 23 such that $\pi^{(p)}\varphi_1 = \pi^{(p)}\varphi_2 = 0$. We then have: $H_0^{(p)}\varphi = 0$ and consequently: $P_t H_0^{(p)}\varphi = 0$ for all $t > 0$. Using $P_t H_0^{(p)}\varphi = H_0^{(p)}P_t\varphi$ (see Davis 1993, e.g.) and Eq. 3, we get $P_t\varphi = \varphi$ and:

$$ \varphi(i, x) = \lim_{t \to +\infty} \frac{1}{t} \int_0^t (P_t \varphi)(i, x) = \pi^{(p)} \varphi $$

for all $(i, x) \in E \times V$. Whence the result, using $\pi^{(p)}\varphi = \pi^{(p)}\varphi_1 - \pi^{(p)}\varphi_2 = 0$. □
*Remark 4* Just as for the transitory results, in the case when the data depend on some family of parameters *P*, the evaluation of all $\lim_{t \to +\infty} \frac{1}{t} \frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t)$ only requires one single computation of both $\pi^{(P)}$ and $Uh^{(P)}$, which seems cheaper than finite differences methods in case of a large *P*.

# 6 Numerical Procedure

In this section, we propose some method for the numerical evaluation of $\frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t)$ and of its asymptotic rate per unit time when $t \to +\infty$, based on the implicit finite volume scheme given in Eymard et al. (2008). These methods are obtained by translating the continuous procedure described above into the discrete setting of finite volume methods, and are provided without proof.

## 6.1 Principle of the Numerical Scheme in the Transient Case

Here, we consider the general case where any quantity (except $\rho_0$) may depend on $p$. We consider a mesh (or partition) $\mathcal{M}$ of $V \subset \mathbb{R}^d$ which satisfies some regularity hypotheses (details are given in Eymard et al. 2008), and some time step $\delta > 0$. Realizing a discrete version of the Chapman-Kolmogorov Eq. 4 and starting from $\rho^{(0)} = \rho_0$, the scheme resumes to compute the family of real values $(\rho_{i,K}^{(n)})_{(i,K,n) \in E \times \mathcal{M} \times \mathbb{N}}$, such that $\rho_{i,K}^{(n)}$ is an approximation of $\int_K \rho_t^{(p)}(i, dx)$ for $t \in [n\delta, (n+1)\delta]$.

Introducing an infinite matrix $A = (A_{(i,K)(j,L)})_{(i,K),(j,L) \in E \times \mathcal{M}}$, the scheme can be written as

$$ \sum_{(j,L) \in E \times \mathcal{M}} A_{(i,K)(j,L)} \rho_{j,L}^{(n+1)} + \frac{\rho_{i,K}^{(n+1)} - \rho_{i,K}^{(n)}}{\delta t} = 0 $$

for all $(i, K, n) \in E \times \mathcal{M} \times \mathbb{N}$ (see the quoted reference for the detailed coefficients), or equivalently

$$ A \rho^{(n+1)} + \frac{\rho^{(n+1)} - \rho^{(n)}}{\delta t} = 0 \qquad (24) $$

for all $n \in \mathbb{N}$, setting $\rho^{(n)} := (\rho_{i,K}^{(n)})_{(i,K) \in E \times \mathcal{M}}$.

For $t = M\delta$ with $M \ge 1$, a discrete approximation of $R_{\rho_0}^{(p)}(t)$ then is $\tilde{R}_M^{(p)}$ with:

$$ \tilde{R}_M^{(p)} = \sum_{n=1}^{M} \sum_{i \in E} \sum_{K \in \mathcal{M}} \& \rho_{i,K}^{(n)} \tilde{h}_{i,K} \qquad (25) $$

and

$$ \tilde{h}_{i,K} = \frac{1}{m(K)} \int_K h^{(p)}(i, x) \, dx \qquad (26) $$

for all $(i, K) \in E \times \mathcal{M}$, where $m(K)$ stands for the $d$-dimensional Lebesgue measure of $K$.
In order to mimic the procedure used to evaluate $\frac{\partial}{\partial p} R_{\rho_0}^{(p)}(t)$ in Eq. 15, we define a discrete version $\tilde{H}$ of $H$:

$$ (\tilde{H}(\theta))^{(n)} = A^t \theta^{(n)} + \frac{\theta^{(n)} - \theta^{(n+1)}}{\delta} \quad (27) $$

for all bounded families $\theta = (\theta_{i,K}^{(n)})_{(i,K,n) \in E \times \mathcal{M} \times \mathbb{N}}$.

We next introduce a discrete approximation of the importance functions $\varphi_t$ as $\bar{\varphi} = (\bar{\varphi}_{i,K}^{(n)})_{(i,K,n) \in E \times \mathcal{M} \times \mathbb{N}}$, single bounded family solution of:

$$ (\bar{H}(\bar{\varphi}))^{(n)} = \bar{h} \quad \text{for all } n \in \{0, \dots, M-1\} \quad (28) $$

with $\bar{\varphi}_{i,K}^{(n)} = 0$ for all $(i, K) \in E \times \mathcal{M}$ and all $n \ge M$. The mathematical study of the well-posedness of Eq. 27 and of the resolution of Eq. 28 might be driven under classical hypothesis.

Following the same calculation steps as in the continuous case, we now get:

$$ \frac{\partial}{\partial p} (\bar{R}_M^{(p)}) = \delta \sum_{i \in E} \sum_{K \in \mathcal{M}} \sum_{n=1}^{M} \rho_{i,K}^{(n)} \left( \frac{\partial}{\partial p} \bar{h}_{i,K} - \left( \frac{\partial \bar{H}}{\partial p} (\bar{\varphi}) \right)_{i,K}^{(n-1)} \right) \quad (29) $$

where, for all $(i, K, n) \in E \times \mathcal{M} \times \mathbb{N}$:

$$ \left( \frac{\partial \bar{H}}{\partial p} (\bar{\varphi}) \right)_{i,K}^{(n)} = \sum_{(j,L) \in E \times \mathcal{M}} \frac{\partial A_{(j,L)(i,K)}}{\partial p} \bar{\varphi}_{j,L}^{(n)}. $$

In the case when the data depends on a family $P$ of parameters, the numerical computation of $\frac{\partial}{\partial p} (R_M^{(P)})$ for all $p \in P$ hence requires to solve two implicit volume schemes provided by Eqs. 24 and 28, which are independent on the choice of $p$ in $P$ (see Remark 3). Due to the definition of $\tilde{H}$ (see Eq. 27), such schemes appear as dual schemes, which clearly simplifies their implementation. Also, the different summations required for each $p \in P$ for the estimation of Eq. 29 are made simultaneously by solving the schemes, which helps saving CPU time and memory stocking size too.

## 6.2 Asymptotical Numerical Procedure

Following the discretization technique introduced in the transient case, a discrete solution of the asymptotic problem is a family of real values $(\bar{\pi}_{i,K})_{(i,K) \in E \times \mathcal{M}}$, such that $\bar{\pi}_{i,K}$ is an approximation of the quantity $\int_K \pi^{(p)}(i, dx)$. Using again the infinite matrix $A = (A_{(i,K)(j,L)})_{(i,K),(j,L) \in E \times \mathcal{M}}$, the asymptotic scheme can then be written as:

$$ \sum_{(j,L) \in E \times \mathcal{M}} A_{(i,K)(j,L)} \bar{\pi}_{j,L} = 0 \quad (30) $$

for all $(i, K) \in E \times \mathcal{M}$, under the constraint

$$ \sum_{(j,L) \in E \times \mathcal{M}} \bar{\pi}_{j,L} = 1. \quad (31) $$
Equivalently, the scheme writes:

$$A\bar{\pi} = 0 \text{ and } \bar{\pi} \cdot \mathbf{1} = 1, \qquad (32)$$

where **1** is the constant family with the required dimension and generic term equal to 1, and where $\cdot$ stands for the dot product $(x \cdot y = \sum_{(i,K) \in E \times \mathcal{M}} x_{i,K} y_{i,K}$ for all $x = (x_{i,K})_{(i,K) \in E \times \mathcal{M}}$ and $y = (y_{i,K})_{(i,K) \in E \times \mathcal{M}}$). This scheme is assumed to have a unique solution.

An approximation $\bar{R}_{\bar{\pi}}^{(p)}$ of $\lim_{t \to +\infty} \frac{1}{t} R_{\rho_0}^{(p)}(t)$ now is:

$$\bar{R}_{\bar{\pi}}^{(p)} = \sum_{i \in E} \sum_{K \in \mathcal{M}} \bar{\pi}_{i,K} \bar{h}_{i,K} = \bar{\pi} \cdot \bar{h},$$

where $\bar{h}_{i,K}$ is provided by Eq. 26 and $\bar{h} = (\bar{h}_{i,K})_{(i,K) \in E \times \mathcal{M}}$.

In the same way as in the transient case, we mimic the continuous procedure and we consider the discrete version $\bar{H}_0$ of $H_0$, given by:

$$(\bar{H}_0(\theta)) = A^t \theta$$

for all bounded families $\theta = (\theta_{i,K})_{(i,K) \in E \times \mathcal{M}}$.

We next introduce a discrete version $\bar{U}$ of the potential function $Uh^{(p)}$, solution to

$$\bar{H}_0 \bar{U} = (\bar{\pi} \cdot \bar{h}) \mathbf{1} - \bar{h} \text{ and } \bar{\pi} \cdot \bar{U} = 0, \qquad (33)$$

which is assumed to have a unique solution. We finally derive the discrete approximation $\frac{\partial \bar{R}_{\bar{\pi}}^{(p)}}{\partial p}$ of $\lim_{t \to +\infty} \frac{1}{t} \frac{\partial R_{\rho_0}^{(p)}}{\partial p}(t)$, which is given by:

$$\frac{\partial \bar{R}_{\bar{\pi}}^{(p)}}{\partial p} = \bar{\pi} \cdot \frac{\partial \bar{h}}{\partial p} + \bar{\pi} \cdot \left( \frac{\partial \bar{H}_0}{\partial p} \bar{U} \right)$$

with

$$\left( \frac{\partial \bar{H}_0}{\partial p} \bar{U} \right)_{i,K} = \sum_{(j,L) \in E \times \mathcal{M}} \frac{\partial A_{(j,L)(i,K)}}{\partial p} \bar{U}_{j,L}$$

for all $(i, K) \in E \times \mathcal{M}$.

The same remarks as in the transitory case are still valid here and in case when the data depends on a family $P$ of parameters, the numerical computation of all $\frac{\partial \bar{R}_{\bar{\pi}}^{(p)}}{\partial p}$ for each $p \in P$ requires to solve exactly two dual volume schemes, provided by Eqs. 32 and 33.

# 7 A First Example

## 7.1 Presentation—Theoretical Results

A single component is considered, which is perfectly and instantaneously repaired at each failure. The distribution of the life duration of the component ($T_1$) is absolutely continuous with respect of Lebesgue measure, with $E(T_1) > 0$. The successive life
durations make a renewal process. The time evolution of the component is described by the process $(X_t)_{t \ge 0}$ where $X_t$ stands for the time elapsed at time $t$ since the last instantaneous repair (the backward recurrence time). There is one single discrete state so that component $I_t$ is here unnecessary. The failure rate for the component at time $t$ is $\lambda(X_t)$ where $\lambda(\cdot)$ is some non negative continuous and bounded function. The process $(X_t)_{t \ge 0}$ is “renewed” after each repair so that $\mu_{(x)}(\text{dy}) = \delta_0(\text{dy})$ and the evolution of $(X_t)_{t \ge 0}$ between renewals is given by $g(x, t) = x + t$.

We are interested in the rate of renewals on $[0, t]$, namely in the quantity $Q(t)$ such that:

$$Q(t) = \frac{R(t)}{t} = \frac{1}{t} \mathbb{E}_0 \left( \int_0^t \lambda(X_s) \, ds \right)$$

where $R(t)$ is the renewal function associated to the underlying renewal process.

The function $\lambda(x)$ is assumed to depend on some parameter $p > 0$ and to meet with $H'_2 (= H_2)$ requirement. (Here, both $\mu_{(x)}(\text{dy})$ and $\mathbf{v}(x)$ are independent of $p$).

Assuming $\mathbb{E}(T_1^{(p)}) < +\infty$, the process is known to have a unique stationary distribution $\pi^{(p)}$ which has the following probability density function (p.d.f.):

$$f_{\pi}^{(p)}(x) = \frac{\mathbb{P}(T_1^{(p)} > x)}{\mathbb{E}(T_1^{(p)})} = \frac{e^{-\int_0^x \lambda^{(p)}(u)du}}{\mathbb{E}(T_1^{(p)})}. \quad (34)$$

Besides:

$$Q^{(p)}(\infty) = \frac{1}{\mathbb{E}_0(T_1^{(p)})} = \frac{1}{\int_0^{+\infty} e^{-\int_0^u \lambda^{(p)}(u)du} du}, \quad (35)$$

which is a direct consequence from the key renewal theorem.

We now provide conditions which ensure the process to be uniformly ergodic ($H_3$):

**Proposition 4** Let us assume that $\mathbb{E}(e^{\delta T_1}) < +\infty$ for some $0 < \delta < 1$ and that $T_1$ is new better than used (NBU), namely such that for all $x, t \ge 0$, we have $\bar{F}(x+t) \le \bar{F}(x)\bar{F}(t)$, where $\bar{F}$ is the survival function $\bar{F}(t) = \mathbb{P}(T_1 > t)$. Then, there are some $C < +\infty$ and $0 < \rho < 1$ such that:

$$|P_t^{(p)} \lambda^{(p)}(x) - \pi^{(p)} \lambda^{(p)}| \le C \rho^t$$

for all $x \in \mathbb{R}_+$.

*Proof* We use the following result by Konstantopoulos and Last (1999): assume that $\mathbb{E}(e^{\delta T_1}) < +\infty$ for some $\delta > 0$. Let $\eta > 0$ be such that $\eta < \delta$ or $\eta \le \delta \wedge 1$ and define

$$W_{\delta, \eta}(x) = 1 + e^{(\delta - \eta)x} \int_{0}^{+\infty} e^{\delta t} \frac{\bar{F}(x+t)}{\bar{F}(x)} dt \ge 1.$$

Then there are positive constants $C_1 < +\infty$ and $\rho < 1$ such that

$$\| P_t^{(p)}(x, \cdot) - \pi^{(p)} \|_{W_{\delta, \eta}} \le C_1 W_{\delta, \eta}(x) \rho^t$$
where

$$ \|v\|_{W_{\delta, \eta}} = \sup_{|g| \le W_{\delta, \eta}} |vg|. $$

Now, under the assumption $\mathbb{E}(e^{\delta T_1}) < +\infty$ for some $0 < \delta < 1$, let us choose $\eta = \delta$. As $T_1$ is NBU, we know that $\bar{F}(x+t) \le \bar{F}(t)$ for all $x, t \ge 0$. Setting $f$ to be the p.d.f. of $T_1$, we derive:

$$ \begin{aligned} W_{\delta, \eta}(x) &\le 1 + \int_0^{+\infty} e^{\delta t} \bar{F}(t) \, dt \\ &= 1 + \int_0^{+\infty} e^{\delta t} \int_t^{+\infty} f(u) \, du \, dt \\ &= 1 + \frac{1}{\delta} (\mathbb{E}(e^{\delta T_1}) - 1) \\ &< +\infty \end{aligned} $$

using Fubini's theorem. From the quoted result, we now derive existence of $C_2 < +\infty$ and $\rho < 1$ such that

$$ \| P_t^{(p)}(x, \cdot) - \pi \|_{W_{\delta, \eta}} \le C_2 \rho^t $$

which easily provides the result due to $|\frac{\lambda^{(p)}}{\|\lambda^{(p)}\|_{\infty}}| \le 1 \le W_{\delta, \eta}$. $\square$

## 7.2 Numerical Results

We assume that $T_1$ is distributed according to some Weibull distribution, which is slightly modified to meet with our assumptions:

$$ \lambda^{(\alpha, \beta)}(x) = \begin{cases} \alpha\beta x^{\beta-1} & \text{if } x < x_0 \\ P_{\alpha, \beta, x_0}(x) & \text{if } x_0 \le x < x_0+2 \\ \alpha\beta (x_0+1)^{\beta-1} & \text{if } x_0+2 \le x \end{cases} $$

where $(\alpha, \beta) \in O = ]0, +\infty[ \times ]2, +\infty[$, $x_0$ is chosen such that $T_1 > x_0$ is a rare event (i.e. $P_0(T_1 > x_0) = e^{-\alpha x_0^\beta}$ is small) and $P_{\alpha, \beta, x_0}(x)$ is some smoothing function which makes $x \mapsto \lambda^{(\alpha, \beta)}(x)$ continuous and non decreasing on $\mathbb{R}_+$. For such a failure rate, it is easy to check that assumptions $H'_2$ and $H'_3$ are true, using Proposition 4.

Taking $(\alpha, \beta) = (10^{-5}, 4)$ and $x_0 = 100$ (which ensures $P_0(T_1 > x_0) \simeq 5 \times 10^{-435}$), we compute $IF_p(t) = \frac{p}{Q^{(p)}(t)} \frac{\partial Q^{(p)}(t)}{\partial p}$ and the asymptotic value $IF_p(\infty)$ for $p \in \{\alpha, \beta\}$ by the finite volume methods (EMR) from Section 6.

For comparison purpose, we also compute such quantities by finite differences (FD) using:

$$ \frac{\partial Q(t)}{\partial p} \approx \frac{1}{\epsilon} (Q^{(p+\epsilon)}(t) - Q^{(p)}(t)) $$

for some small $\epsilon$. As for the transitory FD quantities, we use the algorithm from Mercier (2007) which provides an estimate for the renewal function $R^{(p)}(t)$ and hence for $Q^{(p)}(t) = \frac{R^{(p)}(t)}{t}$. For the asymptotic FD quantities, we use the exact formula Eq. 35.
<table><caption>**Table 1** $IF_\alpha(\infty)$ and $IF_\beta(\infty)$ by finite differences (FD) and by the present method (EMR)</caption><thead><tr><th></th><th>ε</th><th>IF<sub>α</sub>(∞)</th><th>Relative error<br>between FD<br>and EMR</th><th>IF<sub>β</sub>(∞)</th><th>Relative error<br>between FD<br>and EMR</th></tr></thead><tbody><tr><td rowspan="5">FD</td><td>10<sup>-2</sup></td><td>4.625×10<sup>-3</sup></td><td>9.8 × 10<sup>-1</sup></td><td>2.8242</td><td>1.1 × 10<sup>-3</sup></td></tr><tr><td>10<sup>-4</sup></td><td>8.212×10<sup>-2</sup></td><td>6.7 × 10<sup>-1</sup></td><td>2.8214</td><td>1.2 × 10<sup>-4</sup></td></tr><tr><td>10<sup>-6</sup></td><td>2.411×10<sup>-1</sup></td><td>3.5 × 10<sup>-2</sup></td><td>2.8214</td><td>1.1 × 10<sup>-4</sup></td></tr><tr><td>10<sup>-8</sup></td><td>2.499×10<sup>-1</sup></td><td>3.0 × 10<sup>-4</sup></td><td>2.8214</td><td>1.1 × 10<sup>-4</sup></td></tr><tr><td>10<sup>-10</sup></td><td>2.500×10<sup>-1</sup></td><td>7.5 × 10<sup>-5</sup></td><td>2.8214</td><td>1.1 × 10<sup>-4</sup></td></tr><tr><td>EMR</td><td>—</td><td>2.500×10<sup>-1</sup></td><td>—</td><td>2.8211</td><td>—</td></tr></tbody></table>

The asymptotic results are gathered in Table 1.

The comparison between EMR and FD results for small $ε$ clearly validate the method. The results for $IF_β(∞)$ by FD are good and very stable when choosing different values for $ε$. The approximation for $IF_α(∞)$ by FD however requires smaller $ε$ to provide correct results.

The transitory results are next plotted in Figs. 1 and 2 for $t \in [0, 50]$ and different values of $ε$. We observe that FD provides similar results as EMR for $IF_β(t)$ but requires smaller $ε$ for $IF_α(t)$ to get similar results as EMR, just as in the transitory case.

<figure><img src="image_2.png" /><figcaption>**Fig. 1** $IF_α(t)$ by finite differences (FD) for $ε = 10^{-i}$ with $i \in \{2, ..., 7\}$ and by the present method (EMR)</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** $IF_β(t)$ by finite differences (FD) for $ε = 10^{-i}$ with $i \in \{1, 2, 4\}$ and by the present method (EMR)</figcaption></figure>

Finally, comparing $IF_α(t)$ and $IF_β(t)$ for $t \le \infty$, we may note that, for a Weibull distribution, the shape parameter $β$ is much more influent on the rate of renewals than the scale parameter $α$.

# 8 A Second Example

## 8.1 Presentation—Theoretical Results

The following example is very similar to that from Boxma et al. (2005). The main difference is that we here assume $X_t$ to remain bounded ($X_t \in [0, R]$) whereas $X_t$ takes its values in $\mathbb{R}_+$ in the quoted paper.

A tank is considered, which may be filled in or emptied out using a pump. This pump may be in two different states: “in” (state 0) or “out” (state 1). The level of liquid in the tank goes from 0 up to R. The state of the system “pump-tank” at time t is ($I_t, X_t$) where $I_t$ is the discrete state of the pump ($I_t \in \{0, 1\}$) and $X_t$ is the continuous level in the tank ($X_t \in [0, R]$). The transition rate from state 0 (resp. 1) to state 1 (resp. 0) at time t is $\lambda_0(X_t)$ (resp. $\lambda_1(X_t)$). The speed of variation for the liquid level in state 0 is $v_0(x) = r_0(x)$ with $r_0(x) > 0$ for all $x \in [0, R]$ and $r_0(R) = 0$: the level increases in state 0 and tends towards R. Similarly, the speed in state 1 is $v_1(x) = -r_1(x)$ with $r_1(x) > 0$ for all $x \in ]0, R]$ and $r_1(0) = 0$: the level
of liquid decreases in state 1 and tends towards 0. For $i = 0, 1$, the function $\lambda_i$ (respectively $r_i$) is assumed to be continuous (respectively Lipschitz continuous) and consequently bounded on $[0, R]$. The level in the tank is continuous so that $\mu(i, 1 - i, x)(dy) = \delta_x(dy)$ for $i \in \{0, 1\}$ and $x \in [0, R]$.

In order to ensure uniform ergodicity of the process ($\mathcal{H}_3$), we make the following additional assumptions, where, for $i \in E = \{0, 1\}$ and $x, y \in [0, R]$, $\tau_{x,y}^{(i)}$ is the deterministic time for reaching y following the curve $(t, g(i, x, t))$:

$$ \lambda_1(0) > 0 \text{ and } \lambda_0(R) > 0, \qquad (36) $$

$$ \tau_{x_0, R}^{(0)} = \int_{x_0}^{R} \frac{1}{r_0(u)} du = +\infty \qquad (37) $$

$$ \tau_{x_0, 0}^{(1)} = \int_{0}^{y_0} \frac{1}{r_1(u)} du = +\infty \qquad (38) $$

for all $x_0 \in [0, R]$ and all $y_0 \in [0, R]$. This ensures the first jump time $T_1$ to be finite almost surely:

$$ \mathbb{P}_{(0,x_0)} (T_1 < +\infty) = 1 - e^{-\int_{x_0}^{R} \frac{\lambda_0(u)}{r_0(u)} du} = 1 $$

$$ \mathbb{P}_{(1,y_0)} (T_1 < +\infty) = 1 - e^{-\int_{0}^{y_0} \frac{\lambda_1(u)}{r_1(u)} du} = 1 $$

We get the following result:

**Proposition 5** Under assumptions (36–38), the process $(I_t, X_t)_{t \ge 0}$ is positive Harris recurrent with single invariant distribution $\pi$ given by:

$$ \pi(i, dx) = f_i(x) dx $$

for $i = 0, 1$ and

$$ f_0(x) = \frac{K_\pi}{v_0(x)} e^{-\int_{R/2}^{x} \left( \frac{\lambda_1(u)}{v_1(u)} + \frac{\lambda_0(u)}{v_0(u)} \right) du} = \frac{K_\pi}{r_0(x)} e^{\int_{R/2}^{x} \left( \frac{\lambda_1(u)}{r_1(u)} - \frac{\lambda_0(u)}{r_0(u)} \right) du} \qquad (39) $$

$$ f_1(x) = -\frac{K_\pi}{v_1(x)} e^{-\int_{R/2}^{x} \left( \frac{\lambda_1(u)}{v_1(u)} + \frac{\lambda_0(u)}{v_0(u)} \right) du} = \frac{K_\pi}{r_1(x)} e^{\int_{R/2}^{x} \left( \frac{\lambda_1(u)}{r_1(u)} - \frac{\lambda_0(u)}{r_0(u)} \right) du} \qquad (40) $$

for all $x \in [0, R]$, where $K_\pi > 0$ is a normalization constant.

**Remark 5** Though such results are very similar to some special case from Boxma et al. (2005), we have better give here a quick proof due to a few differences in the results, such as some eventual masses for $\pi$ at the bounds of the interval in the quoted paper.

**Proof** Under our assumptions, one can first prove that the process $(I_t, X_t)_{t \ge 0}$ with values in $F := \{0, 1\} \times [0, R]$ is $\varphi$-irreducible (see Meyn and Tweedie 1993b), with maximal irreducibility measure $\varphi = c_{\{0,1\}} \times l$ where $c_{\{0,1\}}$ is the counting measure on $\{0, 1\}$ and $l$ is the Lebesgue measure on $[0, R]$. Besides, the process $(I_t, X_t)_{t \ge 0}$ is non-evanescent, due to values in a compact set. One can also prove that it is a T-process,
namely such that there is some probability measure $a(dt)$ and some kernel $T$ such that for all $i, j \in E$:

$$ \int_{0}^{+\infty} P_t((i, x), (j, A)) a(dt) \geq T((i, x), (j, A)) $$

for all $x \in [0, R]$ and for all Borel set $A \subset [0, R]$, where $T((i, \cdot), (j, A))$ is lower semi-continuous (l.s.c.) for all $(i, j, A)$ and where $T((i, x), E \times [0, R]) > 0$ for all $(i, x)$. Indeed, let $t_0 > 0$ be fixed. For any Borel set $A \subset [0, R]$ with interior $\bar{A}$, and any $(i, j, x) \in E^2 \times [0, R]$, we set

$$ \begin{aligned} T((i, x), (j, A)) &= \mathbf{1}_{\{i=j\}} \mathbb{P}_{(i,x)} (T_1 > t_0) \mathbf{1}_{\bar{A}} (g(i, x, t_0)) \\ &= \mathbf{1}_{\{i=j\}} c_{(i,x)} (t_0) \mathbf{1}_{\bar{A}} (g(i, x, t_0)) \end{aligned} $$

where

$$ c_{(i,x)}(t_0) = e^{-\int_{0}^{t_0} \lambda_i(g(i,x,s))ds} $$

We then have:

$$ \begin{aligned} P_{t_0}((i, x), (j, A)) &= \mathbb{E}_{(i,x)} (I_{t_0} = j, X_{t_0} \in A) \\ &\geq \mathbf{1}_{\{i=j\}} c_{(i,x)} (t_0) \mathbf{1}_{A} (g(i, x, t_0)) \\ &\geq T((i, x), (j, A)) \end{aligned} $$

where $T((i, \cdot), (j, A))$ is l.s.c. because $\bar{A}$ is an open set, so that $(I_t, X_t)_{t \ge 0}$ is a T-process. Using Theorem 3.3 from Meyn and Tweedie (1993b), the process $(I_t, X_t)_{t \ge 0}$ now is Harris recurrent and admits a unique invariant measure $\pi$ up to some multiplicative constant. Besides, $\pi$ and $\varphi = c_{[0,1]} \times l$ are mutually absolutely continuous (see Down et al. 1995) and there is some positive measurable function $f_i$ such that:

$$ \pi(i, dx) = f_i(x) dx $$

Using the fact that $\pi(\cdot, dx)$ is such that $\pi H_0^{(p)} \varphi = 0$ for all $\varphi$ continuously differentiable on $E \times [0, R]$ with

$$ H_0 \varphi (1, x) = -r_1(x) \varphi' (1, x) - \lambda_1(x) \varphi (1, x) + \lambda_1(x) \varphi (0, x) $$

$$ H_0 \varphi (0, x) = r_0(x) \varphi' (0, x) - \lambda_0(x) \varphi (0, x) + \lambda_0(x) \varphi (1, x) $$

one easily finds that

$$ \lambda_{1-i}(x) f_{1-i}(x) - \lambda_i(x) f_i(x) - \frac{d}{dx} (v_i(x) f_i(x)) = 0 $$

for $i = 0, 1$. Solving this system of o.d.e. provides $(f_0, f_1)$ of the form Eqs. 39–40. Checking that, under our assumptions, $\int_0^R f_i(x) dx < \infty$ for $i = 0, 1$, we derive that $\pi$ is a finite measure which can then be normalized in a single way into a probability measure. Consequently, $(I_t, X_t)_{t \ge 0}$ is a positive Harris recurrent process, which ends the proof. $\square$

We now prove that the process $(I_t, X_t)_{t \ge 0}$ is uniformly ergodic ($\mathcal{H}_3$).
**Proposition 6** Under assumptions (36–38), assumption $\mathcal{H}_3$ is true, namely there is some function $f$ fulfilling Eq. 18 such that:

$$|(P_t h)(i, x) - \pi h| \le f(t)$$

for all $(i, x) \in E \times [0, R]$ and all $t \ge 0$.

*Proof* According to Theorem 8.1 by Meyn and Tweedie (1993a), if $(I_t, X_t)_{t \ge 0}$ is bounded in probability on average (condition 1), if there is some $t_0 \ge 0$ such that the kernel $\sum_{n \in \mathbb{N}} P_{t_0 n}((i, x), \cdot)$ possesses an everywhere non trivial continuous component (condition 2) and if the state space is compact (condition 3), then there is some $\rho \in [0, 1]$ such that:

$$\lim_{t \to +\infty} \sup_{(i,x) \in E \times [0, R]} \rho^{-t} \| P_t((i,x), \cdot) - \pi \|_{TV} = 0 \quad (41)$$

where $\|\cdot\|_{TV}$ stands for the norm in total variation. Now, according to Theorem 3.2 from Meyn and Tweedie (1993a), as $(I_t, X_t)_{t \ge 0}$ is a $\varphi$-irreducible $T$-process which is positive Harris recurrent (see Proposition 5 and its proof), the process $(I_t, X_t)_{t \ge 0}$ is bounded in probability on average and condition 1 is true. The second condition is clear too because setting $T_2$ to be the second jump time, we have for any $t_0 > 0$:

$$\begin{aligned} P_{t_0}((i, x), (j, A)) & \ge \mathbb{P}_{(i,x)} (T_1 \le t_0 < T_2, I_{t_0} = j, X_{t_0} \in A) \\ & = \mathbf{1}_{\{j=1-i\}} \int_0^{t_0} c_{(i,x)}(s) \lambda(i, g(i, x, s)) c_{(j,g(i,x,s))}(t_0 - s) \mathbf{1}_A(g(j, g(i, x, s), t_0 - s)) \, ds \end{aligned}$$

for any $(i, j, x) \in E^2 \times [0, R]$ and for any Borel set $A$ (see Cocozza-Thivent et al. 2006b). The kernel $P_{t_0}((i, x), \cdot)$ then possesses an everywhere non trivial continuous component, and $\sum_{n \in \mathbb{N}} P_{t_0 n}((i, x), \cdot)$ does too. As the state space is compact, the three conditions are verified and we derive that Eq. 41 is true. Consequently, there is some $t_1$ such that, for all $t > t_1$, we have:

$$\sup_{(i,x) \in E \times [0, R]} \rho^{-t} \| P_t((i,x), \cdot) - \pi \|_{TV} \le 1$$

and hence:

$$|(P_t h)(i, x) - \pi h| \le \|h\|_{\infty} \rho^t$$

for all $(i, x) \in E \times [0, R]$ and all $t > t_1$. Setting

$$f(t) = \begin{cases} \sup_{(i,x) \in E \times [0, R]} |(P_t h)(i, x) - \pi h| & \text{ift} \le t_1 \\ \|h\|_{\infty} \rho^t & \text{ift} > t_1 \end{cases}$$

easily provides the result.

□
<table><caption>**Table 2** IF<sub>p</sub><sup>(1)</sup> (∞) by FD and EMR</caption><thead><tr><th>p</th><th>FD</th><th>EMR</th><th>Relative error</th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>-3.59 × 10<sup>-2</sup></td><td>-3.57 × 10<sup>-2</sup></td><td>5, 40 × 10<sup>-3</sup></td></tr><tr><td>α<sub>1</sub></td><td>-4.45 × 10<sup>-2</sup></td><td>-4.43 × 10<sup>-2</sup></td><td>3, 65 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>0</sub></td><td>3.19 × 10<sup>-1</sup></td><td>3.17 × 10<sup>-1</sup></td><td>6, 95 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>1</sub></td><td>2.80 × 10<sup>-1</sup></td><td>2.78 × 10<sup>-1</sup></td><td>7, 19 × 10<sup>-3</sup></td></tr><tr><td>a</td><td>4.98 × 10<sup>-1</sup></td><td>4.98 × 10<sup>-1</sup></td><td>1, 06 × 10<sup>-7</sup></td></tr><tr><td>b</td><td>5.09 × 10<sup>-1</sup></td><td>5.09 × 10<sup>-1</sup></td><td>1, 53 × 10<sup>-7</sup></td></tr></tbody></table>

## 8.2 Quantities of Interest

We are interested in two quantities: first, the proportion of time spent by the level in the tank between two fixed bounds $\frac{R}{2} - a$ and $\frac{R}{2} + b$ with $0 < a, b < \frac{R}{2}$ and we set:

$$ 
\begin{aligned}
Q_1(t) &= \frac{1}{t} \mathbb{E}_{\rho_0} \left( \int_0^t \mathbf{1}_{\{\frac{R}{2}-a \le X_s \le \frac{R}{2}+b\}} \mathrm{d}s \right) \\
&= \frac{1}{t} \sum_{i=0}^{1} \int_0^t \int_{\frac{R}{2}-a}^{\frac{R}{2}+b} \rho_s(i, \mathrm{d}x) \mathrm{d}s \\
&= \frac{1}{t} \int_0^t \rho_s h_1 \mathrm{d}s
\end{aligned} 
\quad (42) $$

with $h_1(i, x) = \mathbf{1}_{[\frac{R}{2}-a, \frac{R}{2}+b]}(x)$.

The second quantity of interest is the mean number of times the pump is turned off, namely turned from state “in” (0) to state “out” (1) by unit time, namely:

$$ 
\begin{aligned}
Q_2(t) &= \frac{1}{t} \mathbb{E}_{\rho_0} \left( \sum_{0 \le s \le t} \mathbf{1}_{\{I_{s-}=0 \text{ and } I_s=1\}} \right) \\
&= \frac{1}{t} \mathbb{E}_{\rho_0} \left( \int_0^t \lambda_0(X_s) \mathbf{1}_{\{I_s=0\}} \mathrm{d}s \right) \\
&= \frac{1}{t} \int_0^t \int_0^R \lambda_0(x) \rho_s(0, \mathrm{d}x) \mathrm{d}s \\
&= \frac{1}{t} \int_0^t \rho_s h_2 \mathrm{d}s
\end{aligned} 
\quad (43) $$

with $h_2(i, x) = \mathbf{1}_{\{i=0\}}\lambda_0(x)$.

<table><caption>**Table 3** IF<sub>p</sub><sup>(2)</sup> (∞) by FD and EMR</caption><thead><tr><th>p</th><th>FD</th><th>EMR</th><th>Relative error</th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>-1.81 × 10<sup>-1</sup></td><td>-1.81 × 10<sup>-1</sup></td><td>1, 67 × 10<sup>-4</sup></td></tr><tr><td>α<sub>1</sub></td><td>-1.71 × 10<sup>-1</sup></td><td>-1.71 × 10<sup>-1</sup></td><td>1, 30 × 10<sup>-4</sup></td></tr><tr><td>ρ<sub>0</sub></td><td>-6.22 × 10<sup>-2</sup></td><td>-6.19 × 10<sup>-2</sup></td><td>5, 21 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>1</sub></td><td>-6.05 × 10<sup>-2</sup></td><td>-6.01 × 10<sup>-2</sup></td><td>5, 58 × 10<sup>-3</sup></td></tr></tbody></table>
<table><caption>**Table 4** $IF_p^{(1)}(t)$ for $t=2$ by FD and EMR</caption><thead><tr><th>p</th><th>FD</th><th>EMR</th><th>Relative error</th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>-8.83 × 10<sup>-2</sup></td><td>-8.82 × 10<sup>-2</sup></td><td>1, 08 × 10<sup>-3</sup></td></tr><tr><td>α<sub>1</sub></td><td>-9.10 × 10<sup>-3</sup></td><td>-9.05 × 10<sup>-3</sup></td><td>5, 29 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>0</sub></td><td>4.89 × 10<sup>-1</sup></td><td>4.85 × 10<sup>-1</sup></td><td>7, 51 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>1</sub></td><td>1.97 × 10<sup>-1</sup></td><td>1.97 × 10<sup>-1</sup></td><td>4, 04 × 10<sup>-3</sup></td></tr><tr><td>a</td><td>2.48 × 10<sup>-1</sup></td><td>2.48 × 10<sup>-1</sup></td><td>4, 89 × 10<sup>-4</sup></td></tr><tr><td>b</td><td>7.11 × 10<sup>-1</sup></td><td>7.11 × 10<sup>-1</sup></td><td>7, 77 × 10<sup>-6</sup></td></tr></tbody></table>

## 8.3 Numerical Example

We assume that the system is initially in the state $(I_0, X_0) = (0, R/2)$. Besides, we take:

$$ 
\begin{aligned}
\lambda_0 (x) &= x^{\alpha_0} ; \quad r_0 (x) = (R - x)^{\rho_0} ; \\
\lambda_1 (x) &= (R - x)^{\alpha_1} ; \quad r_1 (x) = x^{\rho_1}
\end{aligned}
$$

for $x \in [0, R]$ with $\alpha_i > 1$ and $\rho_i > 1$. All conditions for irreducibility are here achieved. Our aim is to compute the importance factors with respect to $p$ for $p \in \{\alpha_0, \alpha_1, r_0, r_1, a, b\}$ both in $Q_1(t)$ and $Q_2(t)$, except for parameters $a$ and $b$ which intervenes only in $Q_1(t)$.

We take the following numerical values:

$$ 
\begin{aligned}
\alpha_0 &= 1.05; \quad \rho_0 = 1.2; \quad \alpha_1 = 1.10; \\
\rho_1 &= 1.1; \quad R = 1; \quad a = 0.2; \quad b = 0.2.
\end{aligned}
$$

Similarly as for the first method, we test our results using finite differences (FD). For FD, the transitory results are computed via the finite volume scheme from Eymard et al. (2008) and the asymptotic results via Eqs. 39–40. The results are here rather stable choosing different values for $\varepsilon$ and they are provided for $\varepsilon = 10^{-2}$ in case $p \in \{\alpha_0, \alpha_1, r_0, r_1\}$ and for $\varepsilon = 10^{-3}$ in case $p \in \{a, b\}$. EMR results are computed using the finite volume methods from Section 6. The asymptotic results are given in Tables 2 and 3, and the transitory ones are given in Tables 4 and 5 for $t = 2$.

The results are very similar by FD and EMR both for the asymptotic and transitory quantities, which here again validate the method. Besides, we may observe that the asymptotic results coincides by both methods, even in the case when the velocity field $v(i, x)$ depends on the parameter (here $\rho_i$), which however does not fit with our technical assumptions from Section 5. Due to that (and to other examples where the same remark is valid), one may conjecture that the results from Section 5 are valid under less restrictive assumptions than those given in that section.

<table><caption>**Table 5** $IF_p^{(2)}(t)$ for $t=2$ by FD and EMR</caption><thead><tr><th>p</th><th>FD</th><th>EMR</th><th>Relative error</th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>-2.06 × 10<sup>-1</sup></td><td>-2.06 × 10<sup>-1</sup></td><td>9, 12 × 10<sup>-4</sup></td></tr><tr><td>α<sub>1</sub></td><td>-6.80 × 10<sup>-2</sup></td><td>-6.79 × 10<sup>-2</sup></td><td>2, 12 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>0</sub></td><td>-1.25 × 10<sup>-1</sup></td><td>-1.24 × 10<sup>-1</sup></td><td>4, 27 × 10<sup>-3</sup></td></tr><tr><td>ρ<sub>1</sub></td><td>-4.11 × 10<sup>-3</sup></td><td>-4.03 × 10<sup>-3</sup></td><td>2, 00 × 10<sup>-2</sup></td></tr></tbody></table>
As for the respective importance of the different parameters in $Q_i(t)$, one may note that the importance factors at $t=2$ of $\alpha_0$ and $\rho_0$ in $Q_i(t)$ ($i=1,2$) are clearly higher than the importance factors of $\alpha_1$ and $\rho_1$ in $Q_i(t)$ ($i=1,2$). This must be due to the fact that the system starts from state 0, so that on $[0,2]$, the system spends more time in state 0 than in state 1. The parameters linked to state 0 hence are more important than the ones linked to state 1. Similarly, the level is increasing in state 0 so that the upper bound $b$ is more important than the lower one $a$.

In long-time run, the importance factors of $\alpha_0$ and $\alpha_1$ in $Q_i(t)$ ($i=1,2$) are comparable, which is conform with intuition. The same remark is valid for $\rho_0$ and $\rho_1$, as well as for $a$ and $b$.

Finally, parameters $\rho_0$ and $\rho_1$ are more important than parameters $\alpha_0$ and $\alpha_1$ in $Q_1(t)$, conversely to what happens in $Q_2(t)$. This seems coherent with the fact that quantity $Q_1(t)$ is linked to the level in the tank, and consequently to its evolution, controlled by $\rho_0$ and $\rho_1$, whereas quantity $Q_2(t)$ is linked to the transition rates, and consequently to $\alpha_0$ and $\alpha_1$.

**Acknowledgements** The authors would like to thank Anne Barros, Christophe Bérenguer, Laurence Dieulle and Antoine Grall from Troyes Technological University (Université de Technologie de Troyes) for having drawn their attention to the present subject.

## Appendix

This appendix gives a result of existence, uniqueness and regularity for the solution of the equation:

$$ \varphi(t, i, z) = \varphi_0(i, z) + \int_{t_0}^{t} \mathcal{F}(s, \varphi(s, \cdot, \cdot))(i, z) \, ds $$

We shall use the notations:

- $I = ]a, b[$ bounded interval of $\mathbb{R}$, $E$ finite set and $W$ open subset of $\mathbb{R}^K$,

- $C_1^b(E, W)$ set of bounded function from $E \times W$ to $\mathbb{R}$, continuously differentiable with respect to the variables in $W$ with bounded partial derivatives,

- if $\Psi \in C_1^b(E, W)$, $\nabla\Psi$ is the vector of partial derivatives with respect to the variables in $W$,

- if $\Psi \in C_1^b(E, W)$, $\|\Psi\|_{C_1^b(E, W)} = \max(\|\Psi\|_{\infty}, \|\nabla\Psi\|_{\infty})$,

- $C_1^b(I, E, W)$ set of bounded functions $\varphi(s, i, z)$ from $I \times E \times W$ to $\mathbb{R}$, continuous in the variables in $I \times W$, such that for all $s$ in $I$, the function $(i, z) \to \varphi(s, i, z)$ belongs to $C_1^b(E, W)$ and such that $\sup_{s \in I} \|\varphi(s, \cdot, \cdot)\|_{C_1^b(E, W)} < +\infty$,

- if $\varphi$ belongs to $C_1^b(I, E, W)$, $\|\varphi\|_{C_1^b(I, E, W)} = \sup_{s \in I} \|\varphi(s, \cdot, \cdot)\|_{C_1^b(E, W)}$,

- $C_2^b(E, W)$ set of bounded function from $E \times W$ to $\mathbb{R}$, two times continuously differentiable with respect to the variables in $W$ with bounded partial derivatives,

- if $\Psi \in C_2^b(E, W)$, $J(\Psi)$ is the matrix of second order derivatives with respect to the variables in $W$,

- if $\Psi \in C_2^b(E, W)$, $\|\Psi\|_{C_2^b(E, W)} = \max(\|\Psi\|_{\infty}, \|\nabla\Psi\|_{\infty}, \|J(\Psi)\|_{\infty})$,

- $C_2^b(I, E, W)$ set of bounded functions $\varphi(s, i, z)$ from $I \times E \times W$ to $\mathbb{R}$, continuous in the variables in $I \times W$, such that for all $s$ in $I$, the function $(i, z) \to \varphi(s, i, z)$ belongs to $C_2^b(E, W)$ and such that $\sup_{s \in I} \|\varphi(s, \cdot, \cdot)\|_{C_2^b(E, W)} < +\infty$,
- if $\varphi$ belongs to $C_2^b(I, E, W)$, $\|\varphi\|_{C_2^b(I, E, W)} = \sup_{s \in I} \|\varphi(s, \cdot, \cdot)\|_{C_2^b(E, W)}$.

**Theorem 4** Let us consider the equation

$$ \varphi(t, i, z) = \varphi_0(i, z) + \int_{t_0}^{t} \mathcal{F}(s, \varphi(s, \cdot, \cdot))(i, z) \, ds $$

and let us suppose that:

- $t_0 \in I, t \in I, i \in E, z \in W$,

- $\varphi_0$ belongs to $C_1^b(E, W)$,

- if $\varphi$ belongs to $C_1^b(I, E, W)$, then the application $(s, i, z) \to \mathcal{F}(s, \varphi(s, \cdot, \cdot))(i, z)$ belongs to $C_1^b(I, E, W)$,

- if $\Psi_1$ and $\Psi_2$ are in $C_1^b(E, W)$ there exists a positive constant $C$ such that for all $s$ in $I$:

$$ \begin{aligned} \|\mathcal{F}(s, \Psi_1) - \mathcal{F}(s, \Psi_2)\|_{\infty} &\le C\|\Psi_1 - \Psi_2\|_{C_1^b(E, W)} \\ \|\nabla \mathcal{F}(s, \Psi_1) - \nabla \mathcal{F}(s, \Psi_2)\|_{\infty} &\le C\|\Psi_1 - \Psi_2\|_{C_1^b(E, W)} \end{aligned} $$

Then there exists a unique solution which belongs to $C_1^b(I, E, W)$ and is continuously differentiable with respect to the first variable. If, in addition, we have

- $\varphi_0$ belongs to $C_2^b(E, W)$,

- if $\varphi$ belongs to $C_2^b(I, E, W)$, then the application $(s, i, z) \to \mathcal{F}(s, \varphi(s, \cdot, \cdot))(i, z)$ belongs to $C_2^b(I, E, W)$,

- if $\Psi_1$ and $\Psi_2$ are in $C_2^b(E, W)$ there exists a positive constant $C$ such that for all $s$ in $I$:

$$ \|J(\mathcal{F}(s, \Psi_1)) - J(\mathcal{F}(s, \Psi_2))\|_{\infty} \le C\|\Psi_1 - \Psi_2\|_{C_2^b(E, W)} $$

then there exists a unique solution which belongs to $C_2^b(I, E, W)$ and is continuously differentiable with respect to the first variable.

**Proof** Let us prove the first part of the theorem. To prove the existence of a solution, starting from an initial function $\varphi^{(0)}$ of $C_1^b(I, E, W)$, we define a sequence $\varphi^{(n)}$ of functions of $C_1^b(I, E, W)$ by:

$$ \varphi^{(n+1)}(t, i, z) = \varphi_0(i, z) + \int_{t_0}^{t} \mathcal{F}(s, \varphi^{(n)}(s, \cdot, \cdot))(i, z) \, ds \quad (44) $$

Thanks to theorem of differentiation under integral, we can prove that these functions are in $C_1^b(I, E, W)$ and that the vector of partial derivatives with respect to variables in $W$ satisfies:

$$ \nabla \varphi^{(n+1)}(t, i, z) = \nabla \varphi_0(i, z) + \int_{t_0}^{t} \nabla \mathcal{F}(s, \varphi^{(n)}(s, \cdot, \cdot))(i, z) \, ds $$

We have

$$ \|\varphi^{(n+1)}(t, \cdot, \cdot) - \varphi^{(n)}(t, \cdot, \cdot)\|_{\infty} \le C \int_{t_0}^{t} \|\varphi^{(n)}(s, \cdot, \cdot) - \varphi^{(n-1)}(s, \cdot, \cdot)\|_{C_1^b(E, W)} \, ds $$
and

$$ \|\nabla \varphi^{(n+1)}(t, \cdot, \cdot) - \nabla \varphi^{(n)}(t, \cdot, \cdot)\|_{\infty} \le C \int_{t_0}^{t} \|\varphi^{(n)}(s, \cdot, \cdot) - \varphi^{(n-1)}(s, \cdot, \cdot)\|_{C_1^b(E, W)} \mathrm{d}s $$

so that

$$ \|\varphi^{(n+1)}(t, \cdot, \cdot) - \varphi^{(n)}(t, \cdot, \cdot)\|_{C_1^b(E, W)} \le C \int_{t_0}^{t} \|\varphi^{(n)}(s, \cdot, \cdot) - \varphi^{(n-1)}(s, \cdot, \cdot)\|_{C_1^b(E, W)} \mathrm{d}s $$

Iterating this inequality, we obtain for all $t$ in $I$:

$$ \|\varphi^{(n)}(t, \cdot, \cdot) - \varphi^{(n-1)}(t, \cdot, \cdot)\|_{C_1^b(E, W)} \le K \frac{C^n (t - t_0)^n}{n!} $$

Hence

$$ \|\varphi^{(n)} - \varphi^{(n-1)}\|_{C_1^b(I, E, W)} \le K \frac{C^n (b - a)^n}{n!} $$

Then the sequence $\varphi^{(n)}$ is a Cauchy sequence in the Banach space $C_1^b(I, E, W)$ and it converges to a function $\varphi$ in $C_1^b(I, E, W)$ as $n$ tends to infinity. Taking the limit in Eq. 44, we obtain:

$$ \varphi(t, i, z) = \varphi_0(i, z) + \int_{t_0}^{t} \mathcal{F}(s, \varphi(s, \cdot, \cdot)(i, z)) \mathrm{d}s $$

In order to prove the uniqueness of the solution, let $\varphi_1, \varphi_2 \in C_1^b(I, E, W)$ be two different solutions of the equation. We have:

$$ \varphi_1(t, i, z) - \varphi_2(t, i, z) = \int_{t_0}^{t} (\mathcal{F}(s, \varphi_1(s, \cdot, \cdot)) - \mathcal{F}(s, \varphi_2(s, \cdot, \cdot))) \mathrm{d}s $$

and

$$ \nabla \varphi_1(t, i, z) - \nabla \varphi_2(t, i, z) = \int_{t_0}^{t} (\nabla \mathcal{F}(s, \varphi_1(s, \cdot, \cdot)) - \nabla \mathcal{F}(s, \varphi_2(s, \cdot, \cdot))) \mathrm{d}s $$

Then

$$ \|\varphi_1(t, \cdot, \cdot) - \varphi_2(t, \cdot, \cdot)\|_{C_1^b(E, W)} \le C \int_{t_0}^{t} \|\varphi_1(s, \cdot, \cdot) - \varphi_2(s, \cdot, \cdot)\|_{C_1^b(E, W)} \mathrm{d}s $$

and by iteration

$$ \|\varphi_1(t, \cdot, \cdot) - \varphi_2(t, \cdot, \cdot)\|_{C_1^b(E, W)} \le K \frac{C^n (t - t_0)^n}{n!} $$

Hence for all $n$ in $\mathbb{N}$:

$$ \|\varphi_1 - \varphi_2\|_{C_1^b(E, W)} \le K \frac{C^n (b - a)^n}{n!} $$

The uniqueness follows. We can prove the second part of theorem in the same way.

<figure><img src="image_4.png" /></figure>
References

Boxma O, Kaspi H, Kella O, Perry D (2005) On/off storage systems with state-dependent input, output, and switching rates. Probab Eng Inf Sci 19(1):1-14

Cao X-R, Chen H-F (1997) Perturbation realization, potentials, and sensitivity analysis of Markov processes. IEEE Trans Automat Contr 42(10):1382-1393

Cartan H (1967) Calcul différentiel. Hermann, Paris

Chiquet J, Limnios N (2008) A method to compute the transition function of a piecewise deterministic Markov process with application to reliability. Stat Probab Lett 78:1397-1403

Cocozza-Thivent C, Eymard R, Mercier S (2006a) A finite-volume scheme for dynamic reliability models. IMA J Numer Anal 26(3):446-471

Cocozza-Thivent C, Eymard R, Mercier S, Roussignol M (2006b) Characterization of the marginal distributions of Markov processes used in dynamic reliability. J Appl Math Stoch Anal 1:1-18 (Art. ID 92156)

Costa O, Dufour F (1999) Stability of piecewise-deterministic Markov processes. SIAM J Control Optim 37(5):1483-1502

Costa O, Dufour F (2003) On the Poisson equation for piecewise-deterministic Markov processes. SIAM J Control Optim 42(3):985-1001

Davis MHA (1984) Piecewise-deterministic Markov processes: a general class of nondiffusion stochastic models. J R Stat Soc Ser B 46(3):353-388 (with discussion)

Davis MHA (1993) Markov models and optimization, monographs on statistics and applied probability, vol 49. Chapman and Hall, London

Down D, Meyn SP, Tweedie RL (1995) Exponential and uniform ergodicity of Markov processes. Ann Probab 23(4):1671-1691

Eymard R, Mercier S, Prignet A (2008) An implicit finite volume scheme for a scalar hyperbolic problem with measure data related to piecewise deterministic markov processes. J Comput Appl Math 222:293-323 (Available online 1 Nov 2007)

Gandini A (1990) Importance and sensitivity analysis in assessing system reliability. IEEE Trans Reliab 39(1):61-70

Konstantopoulos T, Last G (1999) On the use of Lyapunov function methods in renewal theory. Stoch Process Appl 79(1):165-178

Lions J-L (1968) Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles. Avant propos de P. Lelong. Dunod, Paris

Mercier S (2007) Discrete random bounds for general random variables and applications to reliability. Eur J Oper Res 177(1):378-405

Mercier S, Roussignol M (2008) Sensitivity estimates in dynamic reliability. In: Bedford T, Quigley J, Walls L, Alkali B, Daneshkhah A, Hardman G (eds) Advances in mathematical modeling for reliability. IOS, Amsterdam, pp. 208-216

Meyn SP, Tweedie RL (1993a) Stability of Markovian processes. II. Continuous-time processes and sampled chains. Adv Appl Probab 25(3):487-517

Meyn SP, Tweedie RL (1993b) Stability of Markovian processes. III. Foster-Lyapunov criteria for continuous-time processes. Adv Appl Probab 25(3):518-548