A new class of survival distribution for degradation processes subject to shocks

Abstract

Many systems experience gradual degradation while simultaneously being exposed to a stream of random shocks of varying magnitudes that eventually cause failure when a shock exceeds the residual strength of the system. In this paper, we present a family of stochastic processes, called shock-degradation processes, that describe this failure mechanism. In our failure model, system strength follows a geometric degradation process. The degradation process itself is any Lévy process, that is, any stochastic process with stationary independent increments. The shock stream is a Fréchet stochastic process, a process derived from the Fréchet extreme-value distribution. Finally, the shock-degradation process is a convolution of the Fréchet shock process and any one of the candidate degradation processes. The system fails at the first occasion when a shock takes system strength across a threshold at zero. The paper presents results for Wiener diffusion processes and gamma processes as examples of Lévy degradation processes. The paper develops key statistical properties of the process model and its survival distribution, including several that are important for its practical application. As the failure mechanism is a first hitting time event, applications that require regression structures fall within the domain of threshold regression methodology.

Introduction

Shock processes have been studied extensively to explain such diverse phenomena as device failures, insurance claims, earth quakes and business bankruptcies, to name but a few. Early work considered shocks arriving according to a stochastic counting process, combined with a probability of surviving the incident shock (Esary and Marshall 1973; Block and Savits 1978). The latter probability might increase or decrease as a function of the cumulative number of shocks. Later studies considered systems governed by bivariate stochastic processes consisting of one process that generated inter-arrival times for shocks and the other that generated the magnitudes of the shocks. The system is assumed to fail when the arriving shock exceeds some critical level (Shanthikumar and Sumita 1983; Gut and Hüsler 1999). The current literature continues to propose many extensions and refinements of shock models. In (Mercier and Hai Ha 2017), for example, a bivariate failure model is proposed for which some shocks may be fatal while others are damaging but not fatal.

Stochastic processes as models for system degradation have also been studied in depth. The fields of application and types of systems studied range widely over business, economics, engineering, health, material science, and many more. To cite a few of the numerous publications on the subject that have appeared in the last two decades, we mention that a Wiener diffusion process and a Gaussian process with drift are used in (Doksum and Normand 1995) to model a key biomarker for patients infected with human immunodeficiency virus (HIV). A Wiener diffusion process is a model for degradation of a self-regulating heating cable in (Whitmore and Schenkelberg 1997). Geometric Brownian motion and gamma processes are models of degradation investigated in (Park and Padgett 2005), with case demonstrations given for electrical resistors and metal fatigue. A good review of degradation models for systems reliability was given in (Ye and Xie 2015).

Many types of systems fail when they receive a shock that is sufficient to break them. For example, optical fiber may break under a momentary strain. On the other hand, in modeling human health, a hip may fracture in a fall, and an acute exacerbation may cause the death of a patient with cystic fibrosis (Castilone et al. 2000; Aaron et al. 2015; He et al. 2015). In this article, we are interested in the convolution of two processes that each contribute to failure. First, an underlying degradation process describes the slow weakening or deterioration of the system through time. Second, a shock process describes a random stream of irregular adverse impacts that are superimposed on the degrading system. Failure occurs when the strength of the system is reduced to the breaking point by the combination of degradation and a shock. We propose a general model which combines stochastic processes for both shocks and degradation to create a doubly stochastic process. We model the degradation process as a Lévy stochastic process (a process having stationary independent increments) and the shock stream by a Fréchet stochastic process (a classical extreme-value generating process). We refer to the convolution of these two processes as a shock-degradation process. The system fails in this model at the first moment that the system strength reaches zero. As this failure time is a first hitting time (FHT) of a threshold by the shock-degradation process, real-world applications of the model can use threshold regression methods to take account of regression covariates (Lee and Whitmore 2006; 2010; Lee et al. 2010).

In the following sections we define the component degradation and shock processes, describe their statistical properties, explain how they are integrated into a single composite process, and explore the nature and form of the survival distribution produced by the composite process.

Composite shock-degradation process

Degradation and system strength processes

The system strength process is denoted by {Y(t),t≥0}. Our model assumes that the initial system strength is Y(0)=y₀>0. Strength Y(t) tends to decline over time in physical systems but in other kinds of systems, such as social, economic and biological systems, strength can fluctuate or even increase over time. The measurement units of strength are those of the application context, such as lung capacity in liters, battery power in watt-hours, or financial reserves in millions of dollars. We refer to ‘degradation of strength’ in our general discussion but the possibility of strengthening through time must be contemplated in some contexts.

The degradation process itself is denoted by {W(t),t≥0} with initial value W(0)=0. The system strength process Y(t) and degradation process W(t) are connected as follows:

$$ Y(t)=y_{0}\exp[W(t)]. $$

(1)

Thus, W(t) is the logarithm of the ratio of system strength at time t to initial system strength and therefore is a measure without physical units. We assume {W(t)} is a stochastic process with stationary independent increments. We also assume that the degradation process has a cumulant generating function (c.g.f.) defined on an open interval ${\boldsymbol {\mathcal {Z}}}$. We denote this function for the process level W(1) at time t=1 by

$$ \ln E\{\exp[uW(1)]\}=\kappa(u). $$

(2)

With stationary independent increments, the c.g.f. of W(t) for any t>0 equals tκ(u).

The connection between system strength and degradation defined in (1) implies that system strength changes in increments that are proportional to the residual strength at any moment. As a result, the system strength Y(t) never reaches a point of zero strength but may approach it asymptotically. Consequently, the system cannot fail through degradation alone. As we show shortly, the failure occurs when shocks exceed the residual strength.

We have assumed that the degradation process {W(t),t≥0,W(0)=0} has stationary independent increments and possesses a c.g.f. κ(u). We are therefore restricting our attention to a subclass of Lévy stochastic processes. This subclass is a large one that includes many processes widely encountered in practical applications. The restriction to Lévy processes that have a c.g.f. eliminates Lévy families whose generating distributions do not possess a c.g.f., such as 1/2-stable processes of Brownian hitting times and Cauchy processes. Our subclass retains most of the important Lévy families such as Wiener processes, gamma processes, inverse Gaussian processes, compound Poisson processes, and additive combinations of these processes. Next, we present a couple of important examples to illustrate the range of types covered by our degradation model:

• Wiener diffusion process. Consider a Wiener diffusion process {W(t),t≥0}, with W(0)=0, having mean parameter μ and variance parameter σ²≥0. A Wiener process consists of independent normally distributed increments. Its sample paths are continuous but tend to meander up and down randomly through time. The c.g.f. of a Wiener process has the form:

  $$ \kappa(u)=u\mu+u^{2}\sigma^{2}/2 \quad \text{for }-\infty < u < \infty. $$

  (3)

  When W(t) is a Wiener process, the strength process Y(t)=y₀ exp[W(t)] is a geometric Wiener process.

  In most applications of the Wiener degradation process, the mean parameter μ is negative, indicating that system strength tends to decline. This situation is usual in physical and engineering systems but applications where μ is positive are encountered in social, economic and health systems. In our general mathematical development, we leave the sign of μ unrestricted.

• Gamma process. Consider a gamma process {X(t),t≥0}, with X(0)=0, having scale parameter η>0 and shape parameter ζ>0. A gamma process consists of independent gamma-distributed increments. It has monotonically increasing sample paths that consist of random steps of random size. For our model here, we take W(t)=−X(t). The negative sign assures that the sample paths of the degradation process {W(t)} are monotonically decreasing. The negative-gamma process finds applications in physical and engineering systems that degrade monotonically and undergo irregular wear and tear through time. A gamma model with monotone increasing sample paths can also be accommodated mathematically but is encountered less frequently in real-world applications. The c.g.f. of W(t) as a negative-gamma process has the form:

  $$ \kappa(u)=\zeta \ln\left(\frac{\eta}{\eta+u}\right) \quad \text{for }\eta+u >0. $$

  (4)

  The corresponding strength process has the following geometric negative-gamma form: Y(t)=y₀ exp[W(t)]=y₀ exp[−X(t)].

A special family of strength processes that is important later as a limiting case is one that has deterministic exponential trajectories, as follows:

• Deterministic exponential process. In this family, system strength follows a deterministic exponential time path of the form Y(t)=y₀ exp(λt) where λ denotes the degradation rate parameter. The c.g.f. for this deterministic process has the degenerate form κ(u)=λu. Note that this exponential function describes actual time decay if λ<0, no decay if λ=0 and system strengthening with time if λ>0.

This deterministic degradation process is obtained as a special instance of the Wiener model if we set σ²=0, in which case λ=μ, the Wiener mean parameter. Likewise, the process is also a limiting case of the negative-gamma degradation process if the ratio −ζ/η is fixed at λ and η is increased without limit.

Fréchet shock process

We look at shock streams through a different lens than previous investigators. In particular, we assume that shocks are generated according to a Fréchet process. This process is a family in the class of max-stable processes. For relevant mathematical background and properties, see (de Haan 1978; 1984; Stoev and Taqqu 2005). By adopting this process, we are not modeling the arrival pattern and magnitudes of individual shocks in an explicit fashion. Instead, the time scale is partitioned into intervals rather than viewed as a continuum of time points. The observed process is the sequence of maximum shocks occurring in the consecutive intervals of the partition. The underlying flow of shocks within each time interval remains latent and uncharacterized - only the maximum shock is potentially observable in each interval.

In mathematical terms, let s>0 be any specified time point and n, any natural number. Consider any partition of the time interval (0,s] into n subintervals (0,t₁],…,(t_n−1,t_n], where 0=t₀<t₁<⋯<t_n=s. We assume that the maximum shock V(0,1) experienced in the process during a unit time interval (0,1] has the cumulative distribution function (c.d.f.) P[V(0,1)≤v)]=G(v), where we refer to G(v) as the generating distribution for the shock process. Let V(t_j−1,t_j)>0 be the maximum shock experienced in subinterval (t_j−1,t_j]. We use V_j as a shorthand notation for V(t_j−1,t_j). The sequence of maxima V₁,…,V_n across the partition of (0,s] are mutually independent, with V_j=V(t_j−1,t_j) having the following c.d.f.:

$$ P(V_{j} \leq v)=G(v)^{t_{j}-t_{j-1}}, j=1,2,\ldots,n. $$

(5)

It follows from (5) that the maximum shock over interval (0,s], namely, V(0,s)= max{V₁,…,V_n}, has the same form of c.d.f. as each of the component terms, namely, P[V(0,s)≤v]=G(v)^s.

In our shock model, the generating distribution G(v) is a Fréchet distribution, which has the following c.d.f.:

$$ G(v)=\exp \left[-\left(\alpha/v \right)^{\beta}\right] \quad \text{for \(v>0\), \(\alpha>0\), \(\beta>0\)}, $$

(6)

where α is a scale parameter and β is a shape parameter. We define a Fréchet shock process as any sequence {V₁,…,V_n} for a specified partition, as just described, where the V_j are generated from (5) using the Fréchet distribution in (6) as the generating distribution.

As already noted, the Fréchet distribution in (6) refers to the maximum shock that will be experienced in one unit of time. The measurement unit of a shock is the same as the measurement unit for system strength. Scale parameter α is the e⁻¹ fractile or 37th percentile of the distribution, irrespective of the value of β. Shape parameter β controls the uniformity of shock magnitudes about this fixed percentile. A distribution with a larger value of β has shocks of more uniform intensity and, thus, a smaller probability of extreme shocks in a given interval. The Fréchet distribution only has finite moments of order smaller than β. Figure 1 shows illustrative plots of the probability density function and cumulative distribution function for the Fréchet shock distribution in (6) when parameter α is set to 1 and β takes the three values 0.5, 1 and 1.5. Observe the long right tail of these distributions, which implies that occasional large shocks will occur. The crossing point of the curves in panel (b) represents the 37th percentile of the distribution, which is fixed at α=1 in the figure.

Fig. 1

a Probability density functions and b cumulative distribution functions for the maximum shock in unit time when α is set to 1 and β is set to values of 0.5 (solid curve), 1 (short-dash curve) and 1.5 (long-dash curve). A larger β is associated with a smaller probability of an extreme shock

The Fréchet distribution is one of the classical extreme value distributions. It is the asymptotic distribution of sample maxima generated from any population distribution that is unbounded to the right and only has finite moments of order k>0 or less. Such populations are said to be of the Cauchy type (Gumbel 1953). A Fréchet distribution of order β is the asymptotic distribution for any Cauchy-type distribution with k=β. The role of the Fréchet distribution as an extreme-value distribution gives a rationale for its use as a shock process. We might imagine a large random sample of shocks being drawn from a given Cauchy-type population in each time interval of the process. The maximum values of such samples will be approximately distributed as a Fréchet distribution. As Cauchy-type distribution families form a large and diverse collection (including the Fréchet distribution itself), this rationale suggests that this kind of shock process has potentially wide application.

Survival distributions of the shock-degradation process

We now combine the shock and degradation processes into a single model. We assume that the shock and degradation processes are independent. With this assumption and the mathematical forms we have chosen for the two component processes, we obtain a conjugate pairing that is mathematically tractable and allows us to generate a wide array of useful mathematical results.

The survival function

In this section we derive an expression for the survival function of our shock-degradation process. We denote the system survival time by S and the sample path of the degradation process W(t) over interval (0,s] by ${\boldsymbol {\mathcal {C}}}=\{W(t):0 \leq t \leq s, W(0)=0\}$.

We start the mathematical development by considering a time partition π={0=t₀<t₁<⋯<t_n=s} of time interval (0,s] and denote the norm of this partition by Δ= maxj(t_j−t_j−1). The Fréchet shock process is described by {V₁,V₂,…,V_n}, where V_j is the maximum shock experienced in subinterval (t_j−1,t_j],j=1,…,n, of time partition π. Next, consider the following local extrema of the degradation process for subintervals of time partition π:

$$ \{(W_{Lj},W_{Uj}), j=1,2,\ldots,n\}, $$

(7)

where W_Lj= inf{W(t):t∈(t_j−1,t_j]} and W_Uj= sup{W(t):t∈(t_j−1,t_j]} for j=1,2,…,n. The pair (W_Lj,W_Uj) are the minimum and maximum degradation levels of the system during subinterval (t_j−1,t_j]. The corresponding pairs of extrema for the strength process are (Y_Lj,Y_Uj) for j=1,2,…,n, where Y_Lj=y₀ exp(W_Lj) and Y_Uj=y₀ exp(W_Uj). These pairs of strength extrema provide the following upper and lower bounds on the conditional probability of system survival through subinterval (t_j−1,t_j], given the system has survived to time t_j−1:

$$ P(V_{j} \leq Y_{Lj})\leq P(S>t_{j}|S>t_{j-1}, {\boldsymbol{\mathcal{C}}}) \leq P(V_{j} \leq Y_{Uj}) \quad \text{for time partition }{\boldsymbol{{\pi}}}. $$

(8)

The lower bound in (8) holds because inequality V_j≤Y_Lj is a sufficient condition for survival. This inequality assures that the maximum shock does not exceed system strength during the jth subinterval even if it occurs at the moment of least strength. The upper bound holds because the complementary inequality V_j>Y_Uj assures failure of the system at some moment during the jth subinterval. Thus, the inequality V_j≤Y_Uj is a necessary condition for survival.

The probability bounds in (8) yield the following bound for system survival beyond time interval (0,s]:

$$\begin{array}{@{}rcl@{}} \prod_{j=1}^{n} \left[ P(V_{j} \leq Y_{Lj})\right] \leq P(S>s|{\boldsymbol{\mathcal{C}}}) \leq \prod_{j=1}^{n} \left[P(V_{j} \leq Y_{Uj})\right] \quad \text{ for time partition }{\boldsymbol{{\pi}}}. \end{array} $$

(9)

As $\phantom {\dot {i}\!}P(V_{j} \leq Y_{Lj})=G(Y_{Lj})^{t_{j}-t_{j-1}}$ and $\phantom {\dot {i}\!}P(V_{j} \leq Y_{Uj})=G(Y_{Uj})^{t_{j}-t_{j-1}}$, we can take logarithms of the two bounds in (9) to obtain the following sums:

$$\begin{array}{@{}rcl@{}} & &L_{{\boldsymbol{{\pi}}}}(s)=\sum_{j=1}^{n} \ln P(V_{j} \leq Y_{Lj})=\sum_{j=1}^{n} (t_{j}-t_{j-1})\ln G(Y_{Lj}), \\ & &U_{{\boldsymbol{{\pi}}}}(s)=\sum_{j=1}^{n} \ln P(V_{j} \leq Y_{Uj})=\sum_{j=1}^{n} (t_{j}-t_{j-1})\ln G(Y_{Uj}). \end{array} $$

(10)

Taking the difference U_π(s)−L_π(s) gives:

$$\begin{array}{@{}rcl@{}} U_{{\boldsymbol{{\pi}}}}(s)-L_{{\boldsymbol{{\pi}}}}(s)&=&\sum_{j=1}^{n} (t_{j}-t_{j-1})[\ln G(Y_{Uj})-\ln G(Y_{Lj})] \\ &=&c \sum_{j=1}^{n} (t_{j}-t_{j-1})[\exp(-\beta W_{Lj})-\exp(-\beta W_{Uj})] \end{array} $$

(11)

where c=(α/y₀)^β>0.

As shown in Appendix A.1 Survival function for a Fréchet shock-degradation process, the sum on the righthand side of (11) converges to zero when we consider a sequence of time partitions π having norms Δ that converge to zero. Therefore, upper bound U_π(s) and lower bound L_π(s) converge in the limit to the same definite integral. This integral equals the conditional log-survival function $\ln P(S>s|{\boldsymbol {\mathcal {C}}})$ for interval (0,s], as follows:

$$\begin{array}{@{}rcl@{}} \ln P(S>s|{\boldsymbol{\mathcal{C}}})={\lim}_{\Delta \rightarrow 0}U_{{\boldsymbol{{\pi}}}}(s)={\lim}_{\Delta \rightarrow 0}L_{{\boldsymbol{{\pi}}}}(s)=\int_{0}^{s} \ln G[Y(t)]dt = -cQ(s), \end{array} $$

(12)

where

$$ c=\left(\frac{\alpha}{y_{0}}\right)^{\beta}\quad \text{and} \quad Q(s)=\int_{0}^{s} e^{-\beta W(t)}dt. $$

(13)

The expression $\ln P(S>s|{\boldsymbol {\mathcal {C}}})$ in (12) denotes the logarithm of the conditional probability of no failure in interval (0,s], where ${\boldsymbol {\mathcal {C}}}=\{W(t):0 \leq t \leq s, W(0)=0\}$ is the realized sample path for system degradation. To obtain the survival function itself, which we denote by $\overline {F}(s)$, we take the expectation of $P(S>s|{\boldsymbol {\mathcal {C}}})$ over all possible degradation sample paths ${\boldsymbol {\mathcal {C}}}$ as follows:

$$ \overline{F}(s) = E_{{\boldsymbol{\mathcal{C}}}}\left[P(S>s|{\boldsymbol{\mathcal{C}}})\right] = E_{{\boldsymbol{\mathcal{C}}}}\left\{\exp\left[-cQ(s)\right]\right\}. $$

(14)

We note for future reference that the right-hand expression in (14) is the Laplace transform of the stochastic integral Q(s) defined in (13), evaluated at c.

The reasoning that leads to survival function (14) shows that the survival time S in our shock-degradation model can be stated compactly in stochastic differential notation as follows:

$$ S=\inf\{t: V(t,t+dt)>Y(t)\}. $$

(15)

The notation V(t,t+dt) represents the maximum shock occurring in the differential time increment (t,t+dt]. Expression (15) states that survival time S is the first moment that the shock exceeds system strength. The conditional log-survival probability in (12) and the survival function in (14) represent a two-step probabilistic evaluation of the stochastic differential expression in (15).

Before proceeding further with the mathematical development, it may be helpful to visualize these shock-degradation processes and the associated failure events when they occur. Figure 2 illustrates sample paths for several types of shock-degradation models. The figure shows simulated paths of system strength Y(t) when a Fréchet shock process is superimposed on several alternative degradation processes W(t). Each illustration has the same simulated shock stream, generated with α=1 and β=1.2. All simulated sample paths have an initial system strength of y₀=20. The panels show, respectively, a system with: (a) constant intrinsic strength (no degradation); (b) deterministic exponential degradation (λ=−0.1), (c) Wiener degradation (μ=−0.02, σ=0.1); and (d) gamma degradation (ζ=2, η=20). In all panels except panel (a), the system happens to fail at a moment between points 15 and 20 on the time scale when a shock first reduces the system strength to zero.

Fig. 2

Simulated paths of system strength Y(t) when a Fréchet shock process is superimposed on several alternative degradation processes W(t). Each illustration has the same simulated shock stream, generated with α=1 and β=1.2, and the same initial strength y₀=20. The panels show, respectively,a system with: a constant intrinsic strength (no degradation); b deterministic exponential degradation (λ=−0.1), c Wiener degradation (μ=−0.02, σ=0.1); and d gamma degradation (ζ=2, η=20)

Lower bound for the survival function

A further result of considerable value is a lower bound for the entire class of survival functions. To derive the lower bound, we have from (14) that:

$$ \ln \overline{F}(s) = \ln E_{{\boldsymbol{\mathcal{C}}}}\left[P(S>s|{\boldsymbol{\mathcal{C}}})\right] \geq E_{{\boldsymbol{\mathcal{C}}}}\left[-cQ(s)\right]=-{cE}_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)\right]. $$

(16)

The result follows from Jensen’s inequality because a logarithmic function is concave downward. Evaluating the right-hand expectation in (16) for any degradation process with stationary independent increments gives the following survival function $\overline {F}_{L}(s)$ as a lower bound for all survival functions in the family:

$$\begin{array}{@{}rcl@{}} \ln \overline{F}_{L}(s)&=& -{cE}_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)\right]=-c\int_{0}^{s} E_{{\boldsymbol{\mathcal{C}}}}\left[e^{-\beta W(t)}\right]dt \\ &=& -c\int_{0}^{s} e^{\kappa t}dt=-c\left[\frac{e^{\kappa s}-1}{\kappa}\right]. \end{array} $$

(17)

Here κ is short-hand for κ(−β), the c.g.f. of the degradation process W(t) evaluated at −β. Parameter κ must be defined for this lower bound to exist. Its existence is not always assured. For example, for a gamma degradation process, we see from (4) that β<η is required for κ=κ(−β) to be defined. In other words, the shape parameter of the shock process cannot be too large relative to the scale parameter of the gamma process.

We find that this lower bound is quite tight for shock-degradation processes where the degradation process is lightly to moderately stochastic. Because the lower bound is often a good approximation and has a tractable closed form, we use it extensively later to explore features of the shock-degradation model.

Numerical evaluation of the survival function

We have not derived a general closed-form expression for the survival function $\overline {F}(s)$ for the shock-degradation family of models. The survival function, however, can usually be evaluated to an adequate approximation for practical applications. The key quantities needed for this evaluation are the expected moments $m_{\ell }(s)=E_{{\boldsymbol {\mathcal {C}}}}\left [Q(s)^{\ell }\right ]$ for ℓ=1,2,… of the stochastic integral Q(s). The exact formula for these moments is a new mathematical result. The formula is derived in Appendix A.2 Exact expressions for expected moments of Q(s) and is given by:

$$ m_{\ell}(s)=E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{\ell}\right] =\ell! \sum_{j=1}^{\ell}\frac{e^{\kappa_{j}s}-1}{\prod_{i=0, i \neq j}^{\ell}(\kappa_{j}-\kappa_{i})} $$

(18)

Here κ_r is shorthand for the value of c.g.f. κ(u) at u=−rβ, for r=0,1,2,…. By definition, κ₀=0. Notice that the lower bound for the survival function given in (17) is obtained directly from the first expected moment $m_{1}(s)=E_{{\boldsymbol {\mathcal {C}}}}\left [Q(s)\right ]$, which comes from (18) by setting ℓ=1 and equating κ₁ with κ in (17).

By calculating successive expected moments m_ℓ(s) from (18), the survival function $\overline {F}(s)$ in (14) can be approximated by finite expansions, as we sketch below and describe in more detail in Appendix A.3 Numerical approximations for the survival function. But caution is required because the larger expected moments of integral Q(s) may grow without limit or even be undefined so the evaluations must be done with care and judgment as we will show. Again, the gamma degradation process offers an example of the difficulty. For the gamma case, quantities κ_r are required to evaluate expected moments of order r or larger from (18) and yet these quantities are only defined if rβ<η. Higher-order expected moments of Q(s) for a Wiener diffusion process also grow without limit.

In the following numerical approximations for the survival function $\overline {F}(s)$ in (14), we temporarily suppress the functional dependence of the random quantities on survival time s to simplify our notation. The approximations are presented in terms of expected central moments of Q which are more accurate. We note that exp(−cQ)= exp(−cm₁) exp[−c(Q−m₁)], where $m_{1}=E_{{\boldsymbol {\mathcal {C}}}}(Q)$ is the first expected moment of Q. Now we replace Q by Q_∗=Q−m₁ and obtain approximations that are based on expected central moments $m^{*}_{\ell }=E_{{\boldsymbol {\mathcal {C}}}}(Q_{*}^{\ell })$ that may be derived from (18) using the usual correspondence of uncentered and centered moments. Observe that $m^{*}_{0}=1$ and $m^{*}_{1}=E_{{\boldsymbol {\mathcal {C}}}}(Q_{*})=0$.

We have found that the following two finite expansions for the survival function $\overline {F}(s)$ are of practical value.

• Taylor series expansion. This approximation is based on expanding the expression for the survival function in a kth order Taylor series and then taking its expectation over ${\boldsymbol {\mathcal {C}}}$, as described in Appendix A.3 Numerical approximations for the survival function. This expansion results in the following approximation:

  $$ \overline{F}_{k}(s)=\exp(-{cm}_{1})\left[\sum_{\ell=0}^{k}(-1)^{\ell}\frac{c^{\ell}}{\ell!}m^{*}_{\ell}\right]. $$

  (19)

• Euler product-limit expansion. This approximation is based on the Euler product limit expansion for exp[−cQ] for which details are given in Appendix A.3 Numerical approximations for the survival function. The kth-order approximation is as follows:

  $$ \overline{F}_{k}(s)=\exp(-{cm}_{1})\left[\sum_{\ell=0}^{k} (-1)^{\ell} \frac{k!}{\ell!(k-\ell)!}\left(\frac{c}{k}\right)^{\ell}m^{*}_{\ell} \right]. $$

  (20)

Because $\exp (-{cm}_{1})=E_{{\boldsymbol {\mathcal {C}}}}\left \{\exp \left [-cQ(s)\right ]\right \}$ is a lower bound for survival function $\overline {F}(s)$, it can be seen that the central-moment approximations are basically adjusting the lower bound multiplicatively to improve the approximation. The lower bound itself is the 1st-order approximation corresponding to k=1 in (19) and (20). As we have shown, the expected moments of Q can diverge or become undefined for larger k but we are free to choose k so this condition is avoided and the error of approximation is kept within acceptable limits.

Panels (a) and (b) of Fig. 3 illustrate survival curves for shock-degradation models having Wiener and gamma degradation processes, respectively. The survival curves have the same parameter values as the simulated sample paths displayed in panels (c) and (d) of Fig. 2. Each panel in Fig. 3 shows a curve for the lower bound in (17), an approximate survival curve calculated from the Taylor series expansion in (19) of order k=4 and, finally, a Kaplan-Meier empirical estimate of the survival curve derived from 1000 simulated survival times. All curves are truncated at survival time s=30. The lower bound curve is quite tight in these cases, especially below the median survival time. We do not show survival curves derived from the Euler product-limit expansion in the figure because these are similar to the Taylor series expansion in these cases for the range of survival times shown.

Fig. 3

Survival curves for two shock-degradation processes. Each illustration has the same shock stream parameters α=1 and β=1.2, and the same initial strength y₀=20. The degradation processes are: a a Wiener degradation process (μ=−0.02, σ=0.1) and b a gamma degradation process (ζ=2, η=20). Each panel shows survival curves corresponding to a lower bound, a 4th order approximation to the true curve based on a Taylor series expansion, and an empirical estimate based on a sample of 1000 simulated cases

Hazard and probability density functions

The hazard function and probability density function are companions of the survival function of a shock-degradation model. Approximations for these functions can be derived mathematically from the corresponding approximations for the survival function. Rather than presenting mathematical forms for these approximations, we find it more convenient to calculate them from the survival function itself using numerical approximations to the differential forms $-d\ln \overline {F}(s)/ds$ and $-d\overline {F}(s)/ds$, respectively.

Because the lower-bound survival function (17) is quite accurate when the degradation process has modest variability, it is instructive to present the mathematical forms of the hazard and probability density functions for this lower bound. The hazard function for the lower bound is:

$$ h_{L}(s)=c\exp(\kappa s) $$

(21)

The hazard is increasing, constant or decreasing exponentially according to whether κ is positive, zero, or negative, respectively. Thus, a roughly exponential shape for the empirical hazard function in a particular application would suggest that the shock-degradation model might be appropriate. Substituting s=0 into (21) gives h_L(0)=c=(α/y₀)^β, which shows quite reasonably that the initial risk of failure depends on α/y₀, the ratio of the 37th percentile of shocks to the initial strength of the system, modulated by the shock shape parameter β.

The probability density function of the lower-bound survival function follows immediately from (17) and hazard function (21). In logarithmic form, it is:

$$ \ln f_{L}(s)=\ln h_{L}(s)+\ln \overline{F}_{L}(s) = \ln(c)+\kappa s-c\left[\frac{e^{\kappa s}-1}{\kappa}\right] $$

(22)

The density function has a positive density of c at the origin and steadily declines as s increases if κ≤c. If κ>c, the density function has a single mode at ln(κ/c)/κ.

Special cases of the survival function

Several special cases of the survival function follow immediately from (14).

• Constant deterministic degradation rate. An immediate result of some interest is obtained in the case where system strength follows a deterministic exponential time path of the form Y(t)=y₀ exp(λt). In this case, the log-survival function can be obtained directly from (14) as:

  $$ \ln \overline{F}(s)=-cQ(s)=-c\int_{0}^{s} e^{-\beta \lambda t}dt=c\left[\frac{e^{-\beta\lambda s}-1}{\beta\lambda}\right] $$

  (23)

  This form of the survival function has already found application in two medical contexts (Aaron et al. 2015; He et al. 2015). By comparing (23) and (17), we see that the lower-bound survival function presented previously is mathematically identical to the survival function for this model with a constant deterministic degradation rate if κ=−βλ in the lower bound formula.

• Pure shock process. A special case of (23) is a pure shock process. For some systems, such as ceramic devices, degradation makes an insignificant contribution to the risk of failure in comparison to the impact of shocks. For this case, the initial strength y₀ is a constant and the log-survival function takes the following form:

  $$ \ln \overline{F}(s)=-cs=-(\alpha/y_{0})^{\beta} s $$

  (24)

  The survival function here is a pure exponential distribution. Although a simple model, its failure rate has a variety of stochastic behaviors depending on the parameters of the shock process and the system strength y₀. To take some limiting cases, for very small β (approaching 0), the survival time is approximately exponentially distributed with a unit rate of decay. For very large β (approaching ∞), failure is almost immediate if y₀<α, is exponentially distributed with a unit rate of decay if y₀=α, and is long delayed if y₀>α.

Joint observation of survival and degradation

Estimation of the shock-degradation model in a practical setting will often involve data that represent longitudinal records on a sample of systems. Each data record consists of periodic readings on the underlying degradation process and ends with either censoring of the record or system failure. This kind of application requires new mathematical results, some of which we now present.

Longitudinal readings on the degradation process

Let us consider more closely the issue of longitudinal readings on the degradation process. We discover that when a degradation process with Fréchet shocks is observed over a brief moment of time, the reading is unlikely to be disturbed significantly by the shock process. Put differently, large shocks occur so infrequently in the time continuum that they are rarely encountered when making one or a few momentary observations on the strength process during a fixed time interval. To see the point mathematically, equation (5) gives the following c.d.f. for the maximum shock V observed in n non-overlapping times intervals, each of length Δt>0:

$$ P(V \leq v)=G(v)^{n\Delta t}=\exp \left[-n \Delta t\left(\frac{\alpha}{v}\right)^{\beta}\right]. $$

(25)

For fixed values of α>0, β>0, v>0 and n, this formula gives a probability that approaches 1 as Δt approaches zero. Thus, in this limiting sense, the maximum shocks observed in n observation intervals can be made arbitrarily small by shrinking the widths of the n observation intervals. Expressed in more practical terms, the theoretical result implies that no material shocks are likely to be present in any finite number of observation intervals if each interval is sufficiently short. The lesson we learn is that material shocks in this kind of process occur sparingly and are rarely discovered at random moments of observation.

Markov decomposition of longitudinal data records

Given the preceding result, we now consider the data elements found in a longitudinal record consisting of periodic readings on system degradation. Our line of development follows the method of Markov decomposition proposed in (Lee and Whitmore 2006). We limit our study now to the lower-bound survival function in (17).

The likelihood function for a longitudinal record will be a product of conditionally independent events of two types. The first type is a failure event. In this event, the system has an initial strength, say y₀, and then fails at time s later. The log-likelihood of this event is given by a log-probability density like that found in (22). The second type is a survival event; more precisely, an event in which the system has an initial strength, say y₀, survives beyond a censoring time s (that is, S>s) and has a recorded strength at time s of, say, y>0. Note that y=Y(s)=y₀ exp[W(s)]=y₀ exp(w), where w denotes the amount of degradation corresponding to strength y. In mathematical notation, the likelihood of this second type of event is:

$$ P\left\{S>s, W(s) \in [w,w+dw]\right\}. $$

(26)

This probability requires derivation of a new result. Its derivation for our shock-degradation process follows the previous line of reasoning. Now, however, the condition previously denoted by set ${\boldsymbol {\mathcal {C}}}$ is replaced by the following more restrictive set for the degradation sample path W(t):

$$ {\boldsymbol{\mathcal{C}}}^{*}=\{W(t):0 \leq t \leq s, W(0)=0, W(s)=w\}. $$

(27)

Set ${\boldsymbol {\mathcal {C}}}^{*}$ represents any sample path of the degradation path that starts at the origin W(0)=0 and ends at level w at time s or, equivalently, any strength sample path that starts at level y₀ at time 0 and ends at level y=y₀ exp(w) at time s. The condition describes a sample path that is pinned down at two end points but is otherwise free to vary.

We note that probability P{S>s,W(s)∈[w,w+dw]} in (26) can be factored into the product of a conditional probability and the p.d.f. g(w) for degradation level w at time s, as follows:

$$ P\left\{S>s, W(s) \in [w,w+dw]\right\}=P\left\{S>s|w\right\}g(w)dw $$

(28)

P.d.f. g(w) is known from the specified form of the degradation process. The conditional survival function P{S>s|w}, which we denote by $\overline {F}(s|w)$, is less straightforward and has yet to be derived in a general form for our shock-degradation model. Here we present a more limited mathematical result.

Survival, conditional on a closing degradation level

For a degradation process with modest variability, we know that a lower bound on the survival function is quite tight. For a pinned degradation process defined by condition ${\boldsymbol {\mathcal {C}}}^{*}$, we have:

$$ \ln \overline{F}(s|w) = \ln E_{{\boldsymbol{\mathcal{C}}}^{*}}\left\{\exp\left[-cQ(s)\right]\right\} \geq -{cE}_{{\boldsymbol{\mathcal{C}}}^{*}}\left[Q(s)\right]=\ln \overline{F}_{L}(s|w). $$

(29)

The inequality follows from Jensen’s inequality when the expectation and logarithmic operators are interchanged. Expectation $E_{{\boldsymbol {\mathcal {C}}}^{*}}\left [Q(s)\right ]$ does not have a general closed form for the family of degradation processes with stationary independent increments; the quantity requires an evaluation of the conditional cumulant generating function $E_{{\boldsymbol {\mathcal {C}}}^{*}}\left \{\exp \left [-\beta W(t)\right ]\right \}$.

We have derived a closed form for the right-hand side of (29) in the important case where the degradation process is a Wiener process. The derivation builds on properties of a Brownian bridge process. The mathematical details of the derivation are given in Appendix Derivation of the survival function conditional on a closing degradation level. Here we simply present the main results. The formula for the lower bound of the survival function for a pinned Wiener process is as follows:

$$ \ln \overline{F}_{L}(s|w)=-c\sqrt{\frac{2\pi s}{\beta^{2}\sigma^{2}}}\exp\left(z_{1}^{2}/2\right)\left[\Phi(z_{2})-\Phi(z_{1})\right], $$

(30)

where

$$z_{1}=\frac{w-\beta(\sigma^{2}/2)s}{\sqrt{\sigma^{2} s}}, \quad z_{2}=\frac{w+\beta(\sigma^{2}/2)s}{\sqrt{\sigma^{2} s}},$$

Φ(·) denotes the standard normal c.d.f.. The p.d.f. of the degradation level W(s) at time s in the Wiener case is given by

$$ g(w)=\frac{1}{\sqrt{2\pi\sigma^{2} s}}\exp[-(w-\mu s)^{2}/2\sigma^{2} s]. $$

(31)

Thus, (30) and (31) are the two components we require for evaluating a tight lower bound for the joint probability in (28), namely, the probability of a system surviving beyond a censoring time s and having the degradation level w at that time.

Other application properties

Practical applications necessarily direct attention to other important properties of shock-degradation processes. We now discuss a couple of these properties.

Probability of no eventual failure (cure rate)

Some versions of shock-degradation processes offer a positive probability that the real-world system that is being modelled will never fail. In medical applications, this probability is often referred to as a cure rate. The cure rate is evident in the lower-bound survival function (17) if parameter κ<0. In this case:

$$ P_{L}(S=\infty)=\exp(c/\kappa) \quad \text{if \(\kappa<0\)}. $$

(32)

To illustrate two contrasting instances of (32), a look at the c.g.f. for a geometric Wiener strength process in (3) shows that κ<0 if μ−βσ²/2>0. Thus, the degradation process in this situation must have μ both positive and large enough to override the variance offset βσ²/2. In this situation, even if the Wiener process has no drift (so μ=0), the survival function will act as if there is a degradation trend toward failure. The tilt towards failure depends on parameter β of the shock distribution as well as the variance parameter σ² of the Wiener process. In contrast, for a gamma degradation process with the c.g.f. defined in (4), we have parameter κ>0 and the degradation sample path decays monotonically and eventual failure is certain so P(S=∞)=0.

Power-transformed time scale

Frequently, an application calls for a power transformation of the survival time scale. For example, with transform s=t^ω, an immediate substitution into (17) gives the following survival distribution on time scale t:

$$ \ln \overline{F}_{P}(t) = \ln \overline{F}_{L}(t^{\omega}) =-c\left[\frac{\exp(\kappa t^{\omega})-1}{\kappa}\right] $$

(33)

This family of survival distributions has already appeared in the distribution theory literature. A survival function with this mathematical form is referred to as an XTGG distribution in (Cordeiro et al. 2016), page 14, for the case where κ>0. Also, (Xie et al. 2002) refers to this type of distribution as a ‘modified Weibull extension.’ The latter article presents some properties, estimation methods, and engineering applications of this type of distribution. What is new here is its connection with the lower bound for survival functions of shock-degradation processes.

Concluding discussion

Our derivation of a family of shock-degradation processes is motivated by the ubiquity of real-world applications in which the strength of a system is degrading stochastically through time and is subjected simultaneously to random shocks of varying size. The system fails when a shock arrives that exceeds its remaining strength. The stochastic process involves a superposition of an intrinsic strength process and an independent shock process. We postulate that system strength decays geometrically according to a degradation process that possesses stationary independent increments and, further, that the shock stream follows a Fréchet process. The combination yields a realistic, tractable and flexible family of processes.

We have derived a number of properties for the shock-degradation process, including a series representation of the survival function and a lower bound for the function that is quite tight in many practical settings. The important situation where survival and the strength of a survivor are jointly observed is considered. For the important case of a Wiener degradation process, we have derived an elegant closed form for the probability of this joint event, which is a new and useful finding. Several additional mathematical results are derived that are important, including an expression for the probability of no eventual failure (which may occur if the system happens to be gradually strengthening rather than degrading). An expression is also obtained for the survival function if the time scale is subject to a power transformation. It is then discovered that the lower-bound survival function, with the transformed time scale, is a modified Weibull extension.

The special case of our shock-degradation model in which the strength process is a deterministic exponential function has been applied already in several studies. This version of the model was used in a study of osteoporotic hip fractures (He et al. 2015). The underlying strength process in the model represented skeletal health. The shock process represented external traumas, such as falls and stumbles, which taken together with chronic osteoporosis might trigger a fracture event. Threshold regression was used to associate time to fracture with baseline covariates of study participants. The same model was considered by (Aaron et al. 2015) to study mortality risks of cystic fibrosis patients. The model was estimated using patient data from a national cystic fibrosis registry. In this setting, acute exacerbations of the disease act as shocks. An exacerbation leads to death on the first occasion when its severity exceeds the respiratory strength of the patient. Practical issues around modelling and estimation for these two studies, as well as for a third study that investigated Norwegian divorces using the same model, were reviewed in the survey paper (Lee and Whitmore 2017).

All of these published applications that employed the simple deterministic-exponential version of the model have used asymptotic maximum likelihood methods for parameter estimation and inferences. The data sets available in these applications were conventional censored survival data. There are important potential applications of the general model, however, that will have more complex types of data than plain censored survival data. For example, in some applications, measurements of both degradation and shock magnitudes prior to failure may be available for model estimation. As an illustration, for patients with chronic obstructive pulmonary disease (COPD), measurements of both lung function and severity levels of acute exacerbations may be gathered through time for each patient. Development of inference methods for the general shock-degradation model based on different data scenarios is a pressing need that we anticipate will be addressed in future research.

Although the shock-degradation model is conceptually clear and well-behaved when the degradation process has only moderate variability, its stochastic behavior becomes chaotic and unpredictable when the degradation process is highly variable. In this situation, the right tail of the survival function becomes heavy and numerical calculation of tail probabilities becomes challenging. Our approximating series expansions are sensitive to this instability in the tails when the degradation process is volatile. An extremely volatile degradation process is not expected in designed or fabricated systems (such as technical devices, engineering structures, or economic and social organizations) but natural physical systems and phenomena frequently have chaotic failure processes. Examples of chaotic natural processes with event times that are extremely unpredictable include, for instance, the time for a cumulus cloud to disperse, an iceberg to collide with a stationary marine structure, or a soap bubble to burst in turbulent air.

A field of application that is the scene of intense current research on stochastic-process models is financial engineering. It seems possible (indeed, likely) that research breakthroughs in stochastic financial models may very well further the development of our shock-degradation model. Business, economic and financial systems are characterized by both stochastic degradation and growth scenarios and these systems are also subject to a wide array of random shocks exactly as in our model. The mathematical advances described in (Hackmann 2015; Hackmann and Kuznetsov 2016), for example, are indicative of this rich source of development potential. Important extensions that will come from future research by these and other investigators might include more efficient numerical methods for calculating tail probabilities of the survival function and other probabilistic quantities for shock-degradation processes.

Appendix

A.1 Survival function for a Fréchet shock-degradation process

A.1.1 Definite integral representation of the conditional log-survival function

Theorem 1

The conditional survival function (s.f.) for a specified degradation sample path${\boldsymbol {\mathcal {C}}}=\{W(t):0 \leq t \leq s, W(0)=0\}$ is given by

$$ P(S>s|{\boldsymbol{\mathcal{C}}})=\exp\left[-cQ(s)\right] $$

(34)

where c=(α/y₀)^β>0 and

$$ Q(s)=\int_{0}^{s}e^{-\beta W(t)}dt. $$

(35)

Proof

We pick up the mathematical development from the main text at expression (11) for the difference U_π(s)−L_π(s):

$$\begin{array}{@{}rcl@{}} U_{{\boldsymbol{{\pi}}}}(s)-L_{{\boldsymbol{{\pi}}}}(s)&=&\sum_{j=1}^{n} (t_{j}-t_{j-1})[\ln G(Y_{Uj})-\ln G(Y_{Lj})] \\ &=&c \sum_{j=1}^{n} (t_{j}-t_{j-1})[\exp(-\beta W_{Lj})-\exp(-\beta W_{Uj})]. \end{array} $$

(36)

It is useful to review briefly the nature of sample paths of our degradation processes, which are Lévy processes. Every Lévy process can be represented as a superposition of processes that may include a Wiener process and independent Poisson processes having varying jump magnitudes. A pure Wiener process (including pure drift) is the only Lévy process with a continuous sample path. All others have discontinuous sample paths because of the Poisson jump component. Jumps of any specified size arrive according to a Poisson process with a size-dependent rate; a rate given by the Lévy measure of the process. For any fixed jump size, say ε>0, the process will experience (almost surely) only a finite number of jumps per unit time that exceed ε but may experience an infinite number of jumps per unit time that are all smaller than ε.

We now examine the limiting value of the difference U_π(s)−L_π(s) when we consider a sequence of time partitions π having norms Δ that converge to zero.

• 1. If our degradation process is a pure Wiener process then difference U_π(s)−L_π(s) will converge to 0 as norm Δ tends to zero because a Wiener process has continuous sample paths. Thus, degradation extremes W_Lj and W_Uj approach each other in each time interval (t_j−1,t_j] as the interval shrinks to zero in the limit.

• 2. If the degradation process has a jump component, we know the largest degradation jump experienced in interval (t_j−1,t_j] of the time partition π cannot exceed $W_{U_{j}}-W_{L_{j}}$. Thus, for any given ε>0, we can divide the time intervals of the time partition π into the following two sets:

  $$J_{0}=\{j: W_{U_{j}}-W_{L_{j}}<\epsilon\}, \quad J_{1}=\{j: W_{U_{j}}-W_{L_{j}} \geq \epsilon\}$$

  We can then express the difference U_π(s)−L_π(s) as the sum A_0π(s)+A_1π(s), where:

  $$A_{0{\boldsymbol{{\pi}}}}(s)=c\sum_{J_{0}} (t_{j}-t_{j-1})[\exp(-\beta W_{Lj})-\exp(-\beta W_{Uj})],$$
  $$A_{1{\boldsymbol{{\pi}}}}(s)=c\sum_{J_{1}} (t_{j}-t_{j-1})[\exp(-\beta W_{Lj})-\exp(-\beta W_{Uj})].$$

  Sum A_1π(s) has a finite number of terms and thus will clearly converge to zero as norm Δ tends to zero. On the other hand, sum A_0π(s) may have an infinite number of terms. We note, however, that:

  $$A_{0{\boldsymbol{{\pi}}}}(s)\leq c [1-\exp(-\beta \epsilon)]\sum_{J_{0}} (t_{j}-t_{j-1})\exp(-\beta W_{Lj})$$

  because $W_{U_{j}}-W_{L_{j}}<\epsilon$ and, therefore:

  $$\exp(-\beta W_{Lj})-\exp(-\beta W_{Uj}) \leq [1-\exp(-\beta \epsilon)]\exp(-\beta W_{Lj}) \quad \text{for all \(j\in J_{0}\)}.$$

  The sum $\sum _{J_{0}} (t_{j}-t_{j-1})\exp (-\beta W_{Lj})$ approaches the definite integral Q(s), defined in (35), as norm Δ approaches 0. This result is evident because time intervals in set J₀ have aggregate measure (0,s] in the limit and W_Lj approaches W(t_j−1) as t_j approaches t_j−1 in the limit, for each interval in J₀. Definite integral Q(s) exists because we have required that c.g.f. κ(u) exist for degradation process {W(t)}. Thus, we have shown that A_0π(s)≤c[1− exp(−βε)]Q(s) in the limit. Finally, as threshold ε>0 can be chosen arbitrarily, we can make it as small as we wish and so drive A_0π(s) to 0 in the limit. In conclusion, therefore, we have shown for a degradation process with a jump component that:

  $${\lim}_{\Delta \rightarrow 0} U_{{\boldsymbol{{\pi}}}}(s)-L_{{\boldsymbol{{\pi}}}}(s)={\lim}_{\Delta \rightarrow 0} A_{0{\boldsymbol{{\pi}}}}(s)+A_{1{\boldsymbol{{\pi}}}}(s)=0.$$

With this last step, the proof is complete. □

A.1.2 Representation of the survival function as an infinite series

The conditional s.f. in (34) can be expanded in an infinite alternating series in successive moments of integral Q(s) as follows:

$$ P(S>s|{\boldsymbol{\mathcal{C}}})=\exp\left[-cQ(s)\right]=\sum_{\ell=0}^{\infty}(-1)^{\ell} \frac{c^{\ell}}{\ell!}Q(s)^{\ell}. $$

(37)

The unconditional survival function is obtained by taking the expectation of $P(S>s|{\boldsymbol {\mathcal {C}}})$ over all sample paths ${\boldsymbol {\mathcal {C}}}$, that is, $\overline {F}(s)=E_{{\boldsymbol {\mathcal {C}}}}\left \{\exp \left [-cQ(s)\right ]\right \}$.

If the expectation operator $E_{{\boldsymbol {\mathcal {C}}}}$ is applied to the infinite alternating series in (37) one is concerned with whether the expectation and summation operations can be interchanged; an issue which invokes the Fubini Theorem or, more precisely, the Fubini-Tonelli Theorem. The interchange would be permitted if either of the following quantities is finite:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left\{\exp\left[cQ(s)\right]\right\} \quad \text{or} \quad \sum_{\ell=0}^{\infty}\frac{c^{\ell}}{\ell!} E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{\ell}\right]. $$

(38)

The first quantity in (38) is the expected value of the reciprocal conditional s.f. in (34), that is, $E_{{\boldsymbol {\mathcal {C}}}}\left [1/P(S>s|{\boldsymbol {\mathcal {C}}})\right ]$. Additionally, from Jensen’s inequality, we have:

$$\overline{F}(s)=E_{{\boldsymbol{\mathcal{C}}}}\left\{\exp\left[-cQ(s)\right]\right\}\geq 1/E_{{\boldsymbol{\mathcal{C}}}}\left\{\exp\left[cQ(s)\right]\right\}=1/E_{{\boldsymbol{\mathcal{C}}}}\left[1/P(S>s|{\boldsymbol{\mathcal{C}}})\right].$$

Neither of the quantities in (38) is guaranteed to be finite in our model setting (as our case examples show) so we must proceed to two finite approximations for the survival function $\overline {F}(s)=E_{{\boldsymbol {\mathcal {C}}}}\left \{\exp \left [-cQ(s)\right ]\right \}$. Both approximations involve evaluations of expected moments of Q(s) for which we derive exact expressions in the next section.

A.2 Exact expressions for expected moments of Q(s)

We have the following expansion for the ℓth expected moment of Q(s) for given survival time s:

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{\ell}\right] &=& E_{{\boldsymbol{\mathcal{C}}}} \left[\prod_{j=1}^{\ell} \int_{0}^{s}e^{-\beta W(t_{j})}{dt}_{j}\right] \\ ~&=& \int_{0}^{s}\cdots\int_{0}^{s} E_{{\boldsymbol{\mathcal{C}}}}\left[\prod_{j=1}^{\ell} e^{-\beta W(t_{j})}\right]{dt}_{1}\cdots {dt}_{\ell} \\ ~&=& \ell! \int_{0}^{s} \int_{0}^{t_{\ell}}\cdots\int_{0}^{t_{2}} E_{{\boldsymbol{\mathcal{C}}}}\left[\prod_{j=1}^{\ell} e^{-\beta W(t_{j})}\right]{dt}_{1}\cdots {dt}_{\ell}. \end{array} $$

(39)

In the last step, the multiple integrals are evaluated over the ordered time points t₁≤t₂≤⋯≤t_ℓ. Because the degradation process W(t) has stationary independent increments over any set of ordered time points, the successive increments W(t_j)−W(t_j−1) for j=1,…,ℓ are independent, where t₀=0. The expectation $E_{{\boldsymbol {\mathcal {C}}}}\left [\prod _{j=1}^{\ell } e^{-\beta W(t_{j})}\right ]$ can be evaluated as follows:

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}\left[\prod_{j=1}^{\ell} e^{-\beta W(t_{j})}\right] &=&E_{{\boldsymbol{\mathcal{C}}}}\left\{\prod_{j=1}^{\ell} e^{-(\ell-j+1)\beta [W(t_{j})-W(t_{j-1})]}\right\} \\ &=&\prod_{j=1}^{\ell} E_{{\boldsymbol{\mathcal{C}}}}\left\{e^{-(\ell-j+1)\beta [W(t_{j})-W(t_{j-1})]}\right\} \\ &=&\prod_{j=1}^{\ell} e^{\kappa_{\ell-j+1}(t_{j}-t_{j-1})} \end{array} $$

(40)

where κ_r is shorthand for the value of c.g.f. κ(u) at u=−rβ, for r=0,1,2,…. By definition, κ₀=0. The derivation requires that the quantities κ_r exist. Their existence is assured if $-r\beta \in {\boldsymbol {\mathcal {Z}}}$ for all natural numbers r of interest. Here ${\boldsymbol {\mathcal {Z}}}$ is the open set on which the c.g.f. κ(u) of our degradation process {Y(t)} is defined.

Substituting (40) into (39) allows the multiple integrals to be evaluated explicitly, giving:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{\ell}\right] =\ell! \sum_{j=1}^{\ell}\frac{e^{\kappa_{j}s}-1}{\prod_{i=0, i \neq j}^{\ell}(\kappa_{j}-\kappa_{i})}. $$

(41)

The following are expressions for the first four moments:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)\right]=\frac{e^{\kappa_{1}s}-1}{\kappa_{1}}, $$

(42)

$$ E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{2}\right]=2!\left[\frac{e^{\kappa_{1} s}-1}{\kappa_{1}(\kappa_{1}-\kappa_{2})}+\frac{e^{\kappa_{2} s}-1}{\kappa_{2}(\kappa_{2}-\kappa_{1})}\right], $$

(43)

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{3}\right]&=& 3! \left[\frac{e^{\kappa_{1} s}-1}{\kappa_{1}(\kappa_{1}-\kappa_{2})(\kappa_{1}-\kappa_{3})}\right. \\ &\left. + \frac{e^{\kappa_{2} s}-1}{\kappa_{2}(\kappa_{2}-\kappa_{1})(\kappa_{2}-\kappa_{3})}+ \frac{e^{\kappa_{3} s}-1}{\kappa_{3}(\kappa_{3}-\kappa_{1})(\kappa_{3}-\kappa_{2})}\right], \end{array} $$

(44)

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}\left[Q(s)^{4}\right]&=& 4!\left[\frac{e^{\kappa_{1} s}-1}{\kappa_{1}(\kappa_{1}-\kappa_{2})(\kappa_{1}-\kappa_{3})(\kappa_{1}-\kappa_{4})}\right. \\ &+& \frac{e^{\kappa_{2} s}-1}{\kappa_{2}(\kappa_{2}-\kappa_{1})(\kappa_{2}-\kappa_{3})(\kappa_{2}-\kappa_{4})} + \frac{e^{\kappa_{3} s}-1}{\kappa_{3}(\kappa_{3}-\kappa_{1})(\kappa_{3}-\kappa_{2})(\kappa_{3}-\kappa_{4})} \\ &\left. + \frac{e^{\kappa_{4} s}-1}{\kappa_{4}(\kappa_{4}-\kappa_{1})(\kappa_{4}-\kappa_{2})(\kappa_{4}-\kappa_{3})}\right]. \end{array} $$

(45)

A.3 Numerical approximations for the survival function

In developing numerical approximations for the survival function $\overline {F}(s)=E_{{\boldsymbol {\mathcal {C}}}}\left \{\exp \left [-cQ(s)\right ]\right \}$, we use notation m_ℓ for $E_{{\boldsymbol {\mathcal {C}}}}\left (Q^{\ell }\right)$, the ℓth expected moment of Q. We also suppress the functional dependence of random quantity Q on survival time s to simplify our notation.

Taylor series expansion

The first approximation to $\overline {F}(s)$ is based on expanding the expression for the survival function in a kth-order Taylor series about zero, followed by taking its expectation over ${\boldsymbol {\mathcal {C}}}$:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left[\exp\left(-cQ\right)\right]=\sum_{\ell=0}^{k}(-1)^{\ell} \frac{c^{\ell}}{\ell!}m_{\ell} +E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(T)}(Q)\right]. $$

(46)

The term $R_{k}^{(T)}(Q)$ is the following Lagrange remainder of the Taylor series expansion:

$$ R_{k}^{(T)}(Q)=(-1)^{k+1}\frac{c^{k+1}}{k!}\int_{0}^{Q} (Q-u)^{k} e^{-cu}du. $$

(47)

Euler product-limit expansion

The second approximation to $\overline {F}(s)$ is based on the expectation of the kth-order Euler product-limit expansion for exp[−cQ], as follows:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left[\exp\left(-cQ\right)\right]=E_{{\boldsymbol{\mathcal{C}}}}\left[\left(1-\frac{cQ}{k}\right)^{k}\right] +E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(E)}(Q)\right], $$

(48)

where $R_{k}^{(E)}(Q)$ denotes the remainder term. By construction, this remainder is:

$$ R_{k}^{(E)}(Q)=\exp(-cQ)-\left[\left(1-\frac{cQ}{k}\right)^{k}\right]. $$

(49)

We note that the first term on the right side of equation (48) can be expanded as the following linear combination of expected moments of Q:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left\{\left[1-\frac{cQ}{k}\right]^{k}\right\}=\sum_{\ell=0}^{k} (-1)^{\ell} \frac{k!}{\ell!(k-\ell)!}\left(\frac{c}{k}\right)^{\ell} m_{\ell}. $$

(50)

Comparing the counterpart terms in (46) and (50) shows that the two approximations differ only in their weights for the first k expected moments of Q.

As both the Taylor series expansion and Euler product-limit expansion equal the same survival function $\overline {F}(s)=E_{{\boldsymbol {\mathcal {C}}}}\left [\exp \left (-cQ\right)\right ]$, it follows by equating (46) and (48) that the expected remainder term $E_{{\boldsymbol {\mathcal {C}}}}\left [R_{k}^{(E)}(Q)\right ]$ in (49) and the one for the Taylor series expansion in (47) are connected mathematically as follows, after some algebraic simplification:

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(E)}(Q)\right]&=&\sum_{\ell=2}^{k}(-1)^{\ell} \frac{c^{\ell}}{\ell!}\left[1-\frac{k!}{(k-\ell)!k^{\ell}}\right]m_{\ell} \\ ~&+&E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(T)}(Q)\right] \quad \text{for }k \geq 2. \end{array} $$

(51)

Approximations based on expected central moments

Re-expression of the preceding approximations in terms of expected central moments of Q sharpens the approximations. We note that

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}[\exp(-cQ)]&=&e^{-{cm}_{1}}E_{{\boldsymbol{\mathcal{C}}}}\{\exp[-c(Q-m_{1})]\} \\ ~&=&e^{-{cm}_{1}}E_{{\boldsymbol{\mathcal{C}}}}[\exp(-{cQ}_{*})]. \end{array} $$

(52)

Here m₁ denotes the first expected moment of Q, that is, $m_{1}=E_{{\boldsymbol {\mathcal {C}}}}(Q)$ and Q_∗=Q−m₁ is the centered Q quantity.

Now we employ the preceding Taylor series expansion and Euler product-limit expansion to expand $E_{{\boldsymbol {\mathcal {C}}}}[\exp (-{cQ}_{*})]$ in (52). We use notation $m^{*}_{\ell }$ for $E_{{\boldsymbol {\mathcal {C}}}}\left (Q_{*}^{\ell }\right)$, the ℓth expected central moment of Q. We therefore have:

$$ E_{{\boldsymbol{\mathcal{C}}}}\left[\exp\left(-cQ\right)\right]=e^{-{cm}_{1}}\sum_{\ell=0}^{k}(-1)^{\ell} \frac{c^{\ell}}{\ell!}m^{*}_{\ell}+e^{-{cm}_{1}}E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(T)}(Q_{*})\right] $$

(53)

$$\begin{array}{@{}rcl@{}} E_{{\boldsymbol{\mathcal{C}}}}\left[\exp\left(-cQ\right)\right]&=&e^{-{cm}_{1}}\sum_{\ell=0}^{k} (-1)^{\ell} \frac{k!}{\ell!(k-\ell)!}\left(\frac{c}{k}\right)^{\ell} m^{*}_{\ell} \\ ~&+&e^{-{cm}_{1}}E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(E)}(Q_{*})\right] \end{array} $$

(54)

By definition, we have $m_{1}^{*}=E_{{\boldsymbol {\mathcal {C}}}}(Q_{*})=0$. Also, because $\exp (-{cm}_{1})=\exp \left \{-{cE}_{{\boldsymbol {\mathcal {C}}}}\left [Q(s)\right ]\right \}$ is a lower bound for survival function $\overline {F}(s)$, it can be seen that the central-moment approximations in (53) and (54) are basically adjusting the lower bound multiplicatively to improve the approximation. In fact, the lower bound itself is the starting approximation with k=1 in the expansion for $E_{{\boldsymbol {\mathcal {C}}}}\left [\exp \left (-cQ^{*}\right)\right ]$.

Approximation bounds for the survival function

The remainder in the Taylor series expansion in (53) has the form:

$$ R_{k}^{(T)}(Q_{*})=(-1)^{k+1}\frac{c^{k+1}}{k!}\int_{0}^{Q_{*}} (Q_{*}-u)^{k} e^{-cu}du. $$

(55)

As the quantity Q_∗∈(−m₁,∞), the term e^−cu in (55) is bounded above by $e^{{cm}_{1}}$ for any u>−m₁ and bounded below by 1−cu for any u. Therefore, we have the following bounds on the remainder term:

$$ (-1)^{k+1}E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(T)}(Q_{*})\right] \leq e^{{cm}_{1}}\frac{c^{k+1}}{(k+1)!}m^{*}_{k+1}. $$

(56)

$$\begin{array}{@{}rcl@{}} (-1)^{k+1}E_{{\boldsymbol{\mathcal{C}}}}\left[R_{k}^{(T)}(Q_{*})\right] &\geq& \frac{c^{k+1}}{k!}\int_{0}^{Q_{*}} (Q_{*}-u)^{k} (1-cu)du \\ ~&=&\frac{c^{k+1}}{(k+1)!}m^{*}_{k+1}-\frac{c^{k+2}}{(k+2)!}m^{*}_{k+2}. \end{array} $$

(57)

As k proceeds through odd and even numbers, the bounds in (56) and (57) become alternately upper and lower bounds and then lower and upper bounds for the remainder term. These bounds on the remainder can be tightened but we will not pursue this line of analysis further.

As anticipated earlier when we discussed the conditions for Fubini’s Theorem to apply, the expected moments of Q can diverge or become undefined for large k and therefore the expected central moments are similarly affected. Fortunately, we are free to choose the finite order of the approximation k as we please and it can be selected so the approximation error is minimized.

We next present a numerical illustration of the approximating Taylor and Euler expansions for the survival function of different orders k to give a concrete sense of how they behave in a typical application. Consider the shock-degradation model for which the shock stream is generated with α=1 and β=1.2 and the degradation process is a Wiener process with μ=−0.02 and σ=0.1. A simulated sample path for this case appears in panel (c) of Fig. 2 and the survival function is approximated in panel (a) of Fig. 3 for survival times up to 30. Table 1 gives selected computational results in this case for $\overline {F}(30)$, the survival probability for s=30. The approximations are stable up to order k=5 but then become unstable. The table shows the Taylor and Euler estimates for the survival probability as well as bounds on these estimates at each order calculated from the bounds for the Taylor remainder term presented earlier. The first line of the table (for order k=1) corresponds to the lower bound $\overline {F}_{L}(30)$ for the survival function at s=30. The table shows this value is 0.2556. As the lower and upper bounds computed at each order k apply to the true probability $\overline {F}(30)$, we know that the largest lower bound across all orders k and the smallest upper bound across all orders k must bound this true probability. The bold entries in the table identify these global lower and upper bounds. The bounds show that $0.2913 \leq \overline {F}(30) \leq 0.2989$.

Table 1 Approximations for the survival probability of an illustrative shock-degradation model for survival time s=30

Derivation of the survival function conditional on a closing degradation level

Consider a set of sample paths for system degradation W(t) that are pinned down at W(0)=0 and W(s)=w. In other words, consider the following set ${\boldsymbol {\mathcal {C}}}^{*}$:

$$ {\boldsymbol{\mathcal{C}}}^{*}=\{W(t):0 \leq t \leq s, W(0)=0, W(s)=w\}. $$

(58)

We now derive the following lower-bound survival function for the shock-degradation family, conditional on the closing degradation level W(s)=w:

$$ \ln \overline{F}_{L}(s|w)=-{cE}_{{\boldsymbol{\mathcal{C}}}^{*}}\left[Q(s)\right]. $$

(59)

Here Q(s) is defined as in (13). Expectation $E_{{\boldsymbol {\mathcal {C}}}^{*}}\left [Q(s)\right ]$ requires an evaluation of the conditional cumulant generating function:

$$\kappa_{t}(-\beta|w)=\ln E_{{\boldsymbol{\mathcal{C}}}^{*}}\left\{\exp\left[-\beta W(t)\right]\right\}.$$

Consider the case where the degradation process is a Wiener process {W(t)} with W(0)=0, as defined in the main text. The derivation builds on properties of a Brownian bridge process. We adapt known results for a standard Brownian bridge to a corresponding pinned Wiener process (see, for example, page 359 of (Ross 1996)). We start with the probability density function for the conditional random variable w_t|w, where W(t)=w_t and W(s)=w, for any 0<t<s:

$$ p(w_{t}|w)=\sqrt{\frac{s}{2\pi\sigma^{2}t(s-t)}}\exp\left\{\frac{-s\left[w_{t}-(t/s)w\right]^{2}}{2\sigma^{2}t(s-t)}\right\}. $$

(60)

As (60) is a normal probability density function, it follows that the corresponding conditional c.g.f. is given by:

$$ \kappa_{t}(-\beta|w)=-\beta(t/s)w+\beta^{2}\sigma^{2}t(s-t)/2s. $$

(61)

The lower bound of the survival function for a pinned Wiener process is therefore derived as follows:

$$\begin{array}{@{}rcl@{}} \ln \overline{F}_{L}(s|w)&=&-c\int_{0}^{s} \exp[\kappa(-\beta|w)]dt \\ ~&=&-c\int_{0}^{s}\exp\left[-\beta(t/s)w+\beta^{2}\sigma^{2}t(s-t)/2s\right]dt \end{array} $$

(62)

The integral in (62) can be evaluated explicitly (see, for example, formula (7.4.32), page 303, of (Abramowitz and Stegun 1964)). After some manipulations, we obtain the following final result:

$$ \ln \overline{F}_{L}(s|w)=-c\sqrt{\frac{2\pi s}{\beta^{2}\sigma^{2}}}\exp(z_{1}^{2}/2)\left[\Phi(z_{2})-\Phi(z_{1})\right] $$

(63)

where

$$z_{1}=\frac{w-\beta(\sigma^{2}/2)s}{\sqrt{\sigma^{2} s}}, \quad z_{2}=\frac{w+\beta(\sigma^{2}/2)s}{\sqrt{\sigma^{2} s}},$$

and Φ(·) denotes the standard normal c.d.f..