Simultaneous Con dence Bands for the Estimation of Expected Discounted Warranty Costs for Coherent Systems under Minimal Repair

This paper develops simultaneous con dence bands using computer intensive methods based on resampling, for the expected discounted warranty costs in coherent systems under physical minimal repair, that is, when the system is observed at its components level and only the component that causes the fault is minimally repaired. For this purpose, it is shown that, under the framework of the Martingale processes and the central limit resampling theorem (CLRT) over stochastic processes, the discounted warranty cost processes satisfy the conditions set by the central limit resampling theorem (CLRT). Additionally, a simulation study is performed on the most relevant variables, such that the nite sample features of the proposed bands may be assessed, based on their actual coverage probabilities. The results in the considered scenarios show that resampling-based simultaneous con dence bands have coverage probabilities that are close to the nominal coverage. In particular, the agreement is good when there are 100 systems or more where a large number of resamples are used for the approximation.


Introduction
When a manufacturer puts a new product on the market, it is expected that a warranty program will come along with the product, which could become a great cost if it is a low-quality product (Thomas 2005).Bai & Pham (2006) highlight that warranty analysis is very important for cost modeling, especially for discounted costs (Blischke & Murthy 1994, Chien 2005, Ja, Kulkarni, Mitra & Patankar 2002, Murthy & Djamaludina 2002, Nguyen & Murthy 1984).Discounted warranty cost models consider the product age and provide a proper measurement of the costs involved in the warranty program, given that they can be treated as random cash ows.Therefore, it is possible to model their evolution throughout the lifetime of the product under warranty, as well as to estimate the necessary fund reserve levels to meet future warranty claims.Several aspects regarding discounted warranty costs and their corresponding reserves have been studied by Bai & Pham (2004), Duchesne & Marri (2009), Ja et al. (2002), Jain & Maheshwari (2006), Patankar & Mitra (1995), and Thomas (1989).
In practice, many products consist of several components, that is, the product can be seen as a system.If every component has its own warranty, they can be combined to produce one warranty for the system, in which it is necessary to keep in mind both the system structure and warranty service costs at its component level (Thomas 1989).Several previous papers have considered warranty cost analysis for component systems (Ritchken & Fuh 1986, Chukova & Dimitrov 1996, Hussain Murthy 1994).The denition of the state of the system immediately before failure depends on the information level one has about the system (Aven & Jensen 1999), so that minimal repairs are classied into two types: statistical minimal repair, which implies replacing the full system for another one just as old as the other one would be if it had not failed (Nguyen & Murthy 1984, Aven 1983, Sheu, Grith & Nakagawa 1995, Ja, Kulkarni, Mitra & Patankar 2001), and physical minimal repair, in which the system is supposed to be observed at its component level and, therefore, only the critical component that caused the system to fail, gets minimal repair (Aven & Castro 2008, González & Bueno 2011).González & Bueno (2011) propose a Martingale estimator for the expected discounted warranty costs for a coherent system under physical minimal repair and include the calculation of specic condence limits which do not make up a simultaneous condence band, given that the aforementioned set of limits generally has no correct coverage probability (Fleming & Harrington 1991).The main purpose of constructing simultaneous condence bands is to assess an estimator's precision, which can be described by the distribution (or a function of it) of that estimator's deviations from its real value (Belyaev 2007).The problem is that the aforementioned distribution is unknown, even if asymptotic convergence results of distribution can be obtained (González & Bueno 2011).In the practice, sample sizes are not always large enough for those approximations to work properly.In general, computer intensive (CI) methods provide a way to nd asymptotically precise approximations of the estimator deviation distributions from the real unknown parameters (Belyaev 2000).The bootstrap, introduced by Efron (1979) is a rather universal method, nevertheless, the need to nd a proper estimator of the real parameter which can describe data distribution, may be a dicult problem, which is why resampling can be used alternatively (Belyaev 2000).
Resampling is used in this paper to develop simultaneous condence bands for the mean function of the discounted warranty costs of a system under physical minimal repair.For this, based on the theoretical framework of Martingale process and the Central Limit of Resampling Theorem (CLRT) over stochastic processes proposed by Belyaev (2000) and Belyaev & Seleznjev (2000), proof is presented that the discounted warranty costs processes under the general lifetime model comply with CLRT conditions.In addition, a simulation study was conducted on the most relevant variables to test the nite sample features of the proposed simultaneous condence bands by means of the actual coverage probability.
Section 2 presents the theoretical framework that is necessary in the development of this paper.The proposal of constructing simultaneous condence bands is developed in Section 3. In Section 4, the performance of the proposed simultaneous condence bands is assessed by means of a simulation study.Section 5 presents the paper's most important conclusions and recommendations.

Theoretical Framework
In this paper it is assumed that there is a coherent system under physical minimal repair, that is, under the physical approach of the general failure model Revista Colombiana de Estadística 41 (2018) 130 (Aven & Jensen 1999).Next section summarizes some theoretical results which are necessary for the development of the remaining sections in the paper.

Physical Minimal Repair Model for a Coherent System and Discounted Warranty Costs
Let us suppose a system with m components, where T is the system lifetime, S i is the lifetime of component i, i = 1, . . ., m and Ñt is the number of minimal repairs of the system in the interval [0, t], dened on a complete probability space (Ω, F, P ) with the ltration F = (F t ) t≥0 , where I(A) is the indicator of the event A. Therefore, the system repair/failure process is observed at the level of its m components.Suppose the following conditions are hold: a) 0 < S i < ∞ P -a.s., i = 1, . . ., m, where P -a.s.denotes that an event E happens almost surely with respect to the probability P .b) For every i = j, P (S i = S j ) = 0, that is, there are not two components failing simultaneously.c) All lifetimes S i are totally inaccessible F t -stopping times, and consequently all the compensators A i of the respective simple counting processes N i t = I (S i ≤ t), are continuous P -a.s., with the Doob-Meyer decomposition Under these assumptions and according to Aven & Jensen (1999) can be shown that where Y i is the critical level of component i, that is, the rst instant after which failure of component i causes the system failure.González & Bueno (2011) show that the system failure indicator process N t = I (T ≤ t), given by has the compensator A t as follows where λ i (t) is the failure rate function for component i.
From (4), it is clear that the F t -intensity of the system on {T > t} is given by (5) Revista Colombiana de Estadística 41 (2018) 130 Let N i t be the number of minimal repairs on the component i in the interval [0, t] and let H i (t) be the minimal repair discounted cost of component i at time t, a deterministic, continuous, decreasing, bounded, and integrable function on the interval [0, t], such that ) be the accumulated cost process by minimal repairs on component i on the interval [0, t], where S ij is the time of the j-th minimal repair of component i and S i1 = S i .González & Bueno (2011) show that B i t has and that when only the critical component causing the system failure is minimally repaired, that is, considering the set Denition 1. (González & Bueno 2011) For each ω ∈ Ω, the set of components which survive its critical level is dened as For each i = 1, . . ., m dene C i (ω) = Then, the minimal repair counting process of the coherent system in [0, t] is and the warranty cost process is given by B Then, from Denition 1 the expected cost of minimal repairs carried out over with Based on the results above, next section gives an estimation method for B * (t) based on a random sample of n identical systems (or n independent copies of the process)., C i(j) , i = 1, . . ., m, j = 1, . . ., n .For each j, let C Φ(j) be the set of critical components for the j-th system, which is dened as in (8), then the minimal repair costs for the j-th system is whose compensator process is From the n independent copies, the following mean processes are obtained, t .
The following are results for the mean processes dened above.
ii.A consistent and unbiased estimator for Var B (n) iii.An approximate 100 (1 − γ) % condence interval for B * (t), is where Z γ/2 is the (1 − γ/2) quantile of the standard normal distribution and Despite having a set of approximately 100 (1 − γ) % level pointwise condence limits in [0, t], given by (13), they do not form 100 (1 − γ) % level simultaneous condence bands.The following section introduces a general denition of simultaneous condence bands for random functions.
Except for functions having a very simple structure, there are no simultaneous condence bands with an exact (1 − γ) coverage probability (Knowles 1988).
Let f be an estimator of the function f (t), based on a random sample of size n, then the weak convergence of processes with the form provides a general method for calculating simultaneous condence bands for the function f (t) (Fleming & Harrington 1991).When where q γ (t) is the (1 − γ) quantile in the distribution of sup 0≤s≤t |Q (s)|, then, asymptotically, Then, the construction of simultaneous condence bands is based on nding q γ (t) which satises the desirable coverage probability on the interval [0, t].
Simultaneous condence bands based on sup 0≤s≤t , will be useful only when sucient conditions for the convergence of reasonable intervals [0, t] are not too restrictive, when q γ (t) is easy to calculate, and when the resulting bands have appealing properties (Fleming & Harrington 1991).Even when the general conditions for weak convergence could be fullled, calculating q γ (t) requires determining the limit process sup 0≤s≤t |Q (s)| to which the process sup 0≤s≤t √ n f s − f (s) converges, and this is not easy when only a sample of n systems is available.
The following section presents the weak approach of processes introduced by Belyaev (2000) and Belyaev & Seleznjev (2000) as an extension of weak convergence of processes, which justies the use of resampling in the approximation of asymptotic distributions.

Weakly Approaching Distributions
Let {L (U n )} n≥1 and {L (V n )} n≥1 be two sequences of distributions of random variables U n and V n which have values on H n , a metric space with metric d n , and let C b (H n , d n ) be the set of all bounded real-valued continuous functions on H n .Denition 2. (Belyaev 2000) The distributions {L (U n )} n≥1 weakly approach Let W n be a random variable (in a metric space W n ) dened on the same probability space of U n previously dened.Suppose that the regular conditional Lemma 1. (Belyaev 2000) Let U n , W n and V n be as dened before.Suppose that ←→ L (V n ) and let Z n be an H n -valued random variable dened on the same probability space of U n , such that The notion of weakly approaching establishes a variant of Lyapunov's Central Limit Theorem (CLT) for R k , as follows.Let U n = {U 1n , . . ., U nn }, n = 1, 2, . .., be a triangular scheme of independent vector-valued random variables, where for each n, U in = (U 1in , . . ., U kin ) T .Let C in , i = 1, . . ., n be the covariance matrix between the k variables of U in and dene U Theorem 2. (Belyaev 2000, CLT in R k ) Suppose that for some constants δ > 0 where Now consider a triangular scheme U n = {U 1n , . . ., U nn } of independent vectorvalued random variables with U in ∈ R k , (i, n) ∈ T , where T = {(i , n ) i = 1, . . ., n, n = 1, 2, . ..}.Let J 1n , . . ., J nn be the indexes of a resample from U n , indicating which of the observations U in is chosen as the i-th element in the resample.Let Revista Colombiana de Estadística 41 (2018) 130 N hn = n i=1 I (J in = h) be the number of times that the observation h is drawn in the resample and dene the vector-valued random variable U 0 which can be interpreted as a centered sum of values, n times randomly sampled with replacement from the components in the statistical data U n .
Theorem 3. (Belyaev 2000, CLRT in R k ) Suppose the assumptions of Theorem 2 are fullled and that E [U hin ] = µ hn , that is, the expectation does not dependent of i.Then, it holds that In general, the resampling process consists of simulating G copies {J r 1n , . . ., J r nn }, r = 1, . . ., G, (which are used to approximate L U 0 Belyaev (2000) shows that to assess the accuracy of some nonparametric estimators, it is necessary to consider more general spaces than R k .For instance, for many non-parametric statistical models it is necessary to consider the Skorokhod space D [0, t] , t > 0, of the so-called cadlag functions.That is, the set of all real functions v (•) dened in [0, t] such that for all s ∈ (0, t) there are limit values v (s ± 0) = lim h↓0 v (s ± h) and v (s) = v (s + 0).

Simultaneous Condence Bands for the Expected Discounted Warranty Cost
In constructing simultaneous condence bands for B * (t) it is important to assess the distribution of the unobservable processes 4 and the use of resampling allow approximating relevant distribution when information about an initial sample of n systems is available.For this purpose, dene where, are stochastic processes in D[0, t], which can be arranged in a triangular scheme U n .
To use the CLRT in D[0, t], it must be veried that the processes U jn (t) satisfy the conditions C-1, C-2 and C-3 established in Theorem 4, which is summarized in the following theorem (see proof in the Appendix B) Theorem 5. Let U jn (t) be as dened in (25).Then for every The following corollary formalizes the application of the CLRT to warranty cost processes.
Corollary 1.Consider a sample of n independent copies of and the triangular scheme U n (t) = {U n1 (t) , . . ., U nn (t)} with U jn (t) given in (25).Then, for the process one can show that Revista Colombiana de Estadística 41 (2018) 130 Proof .By applying Theorem 5 over the process U jn (t), the necessary conditions for Theorem 4 are obtained, from whose application the result is obtained.
In practice the processes U jn (t) are unknown and they need to be estimated.The following section uses estimations of the processes U jn (t) to construct simultaneous condence bands based on resampling.

A Proposal for Simultaneous Condence Bands for the Expected Warranty Cost
The unobservable process U jn (t) can be rewritten as follows: where The following result establishes that the weakly approaching of the process U 0 •n (t), given in ( 26), is kept in an estimated version of itself.
Corollary 2. The process U 0 Proof .By using (28) in ( 26), the identity is obtained, where U 0 Notice that U 0 •n (t) ≡ 0 (since U 0 •n (t) does not dependent of j by denition and n j=1 N jn = n).Then, U 0 •n (t) P −→ 0, therefore by using Lemma 1 the result is obtained.
The latter allows proposing simultaneous condence bands for the expected discounted warranty cost in coherent systems under physical minimal repair.
Theorem 6.An approximate 100 (1 − γ) % simultaneous condence band for B * (t), the expected discounted warranty cost process for a coherent system under physical minimal repair on the interval [0, t], is where, q γ (t) is the (1 − γ) quantile of the empirical distribution of sup 0≤s≤t U 0 •n (s) .
Revista Colombiana de Estadística 41 (2018) 130 Proof .González & Bueno (2011) showed that the process U •n (t) weakly converges in the space D [0, t] to a Gaussian stochastic process denoted by Q (t) which implies By following the denition of simultaneous condence bands introduced in Section 2.3 and from (32) the idea is nding q γ (t) such that: where the process sup Thus, instead of nding q γ (t) satisfying ( 33), the purpose is obtaining q γ (t) such that: This is equivalent, using Corollary 2, to nd the value of q γ (t) such that Thus, q γ (t) can be chosen as the (1 − γ) quantile of the empirical distribution based on resampling from sup 0≤s≤t U 0 •n (s) .Thus, considering that where then, with a value of M large enough, it holds that such that B (n) The following are the factors considered in the simulation study:

Simulation Study
• System type.A three-component parallel system and a 2-out-of-4 components system are considered.To avoid confusion, both systems are treated as k-out-of-m systems, which are denoted by: a) Φ 1:3 t for the three-components parallel system or 1-out-of-3 components system, and b) Φ 2:4 t for the 2-outof-4 components system.
• Number of systems under warranty.This is denoted by n and corresponds to the number of independent and identical copies of the repair/cost process used for constructing simultaneous condence bands.The levels considered are n = 10, 30, 50, 100, 500 and 1000.
• Discount function.Denoted by H i (t) describes the consumer share of physical minimal repair costs for the system during the W warranty term.H i (t) = c i e −t , c i 1 − tW −1 e −t , i = 1, . . ., m were used.
• Number of resamples.This is denoted by G.The levels of G are 500, 1000, 5000 and 10000.
• Partition size.This is denoted by M and it determines how thin is the partition of the warranty period, for the approximation of the supreme of the limit process given in (32).The levels of M are 100, 500 and 1000.
Table 1 summarizes the levels considered in the simulation factors.The following are xed values or simulation parameters: • Component failure distributions.For the systems considered (Φ 1:3 t and Φ 2:4 t ) in each component i, whose associated lifetime is denoted S i , the respective cumulative distribution function F i is considered, thus: • Warranty term.W denotes the time period within which the system is under warranty.The simulation uses W = 5.This can be interpreted as representing ve years or ve thousand use cycles.

Parameter
Fixed value The purpose, at this point, is to assess the performance of the simultaneous condence band proposed in each scenario, by estimating the coverage probabilities.

Coverage Probabilities
Let SCB (n) t be a simultaneous condence band for the function B * (t) in [0, t], based on a sample of n systems, then the coverage probability (CP) for SCB (n) t is dened as: The following is the procedure followed during the simulations: i.For each scenario, generate N simulations of n systems under warranty.
iii.For each scenario, calculate the actual coverage probability for the simultaneous condence band, where the indicator variables determine if the function B * (t) is totally contained within the resulting bands in each simulation.
Since the function of the expected cost B * (t) is unknown, it is approximated with B (n) t for a sample of 100000 systems.Figure 1 shows the functions B * (t) for coherent systems and discount functions under study.
a) 1-out-of-3 components system b) 2-out-of-4 components system Then, the actual coverage probability for the proposed simultaneous condence bands is obtained from (39).
The following are the results of the simulation study for each of the coherent systems considered.

Actual Coverage Probabilities for the 1-out-of-3 Components System
For each value considered of the partition size M of the warranty period, Figure 2 shows the results for analyzing the eects of the size of resamples G and the size of sample n, over the actual coverage probabilities for both discount functions., for the 1-out-of-3 components system, varying G for both discount functions.
Note that for the values of M = 100 and 500 considered for the partition size, there are just small dierences between the CP curves for the resample sizes G studied in each discounted cost model considered.This suggests that for this study the resample size G does not aect the behavior of the actual coverage probabilities of the proposed bands.Nevertheless when dealing with ne partitions M = 1000, the actual coverage probabilities increase because the larger number of resamples.
Figure 2 shows that the dierences in actual coverage probabilities achieved by the two discount functions decrease when the number of systems increases, reaching values close to the nominal level of (1 − γ) used in the simulation.4.4.Actual Coverage Probabilities for the 2-out-of-4 Components System Figure 3 shows the eect of the resample size G over the actual coverage probabilities under dierent sample sizes, for both discount functions., for the 1-out-of-3 components system, varying G for both discount functions.
Figure 3 shows that, similar to, the 1-out-of-3 system, when the number of systems under warranty increases, the actual coverage probabilities of the simultaneous condence bands increase toward the level (1 − γ).At each level of the partition size M , the smaller coverage probabilities suggest an improvement when the resample size G increases.It is worth noting that there are large differences between the coverage probabilities for both discount functions when the number of systems is smaller or equal than 100.But similar values to the nominal (1 − γ) level are achieved in both discount functions when the number of systems is greater than 100.

Examples of Simultaneous Condence Bands
Simultaneous condence bands are illustrated with simulated data for both the 1-out-of-3 components system and the 4-out-of-2 components system for xed values: W = 5; α = 0.05; Figures 4 and 5 show the precision changes within the simultaneous condence bands when the number of systems n under warranty varies, for the 1-out-of-3 and 2-out-of-4 components systems, respectively.
Notice that for both systems in the study, the number of systems under warranty is a key factor in the precision of the proposed simultaneous condence bands.Observe that the bands are narrow when n ≥ 500.Few systems under warranty (n ≤ 100) generate high variability of the estimation of the expected discounted warranty cost as reected in wider simultaneous condence bands.

Conclusions and Future Research
For some statistical models, assessment of the precision of the statistical inferences may be carried out by means of intensive computer methods (Efron 1979, Davison & Hinkley 1997, Belyaev 2007).Resampling was eciently used in this work to obtain simultaneous condence bands, for the expected discounted warranty cost under physical minimal repair.This is a useful tool to assess the precision of the estimator, avoiding the complications of the asymptotic analysis of the related stochastic processes.For instance, it avoids nding the specic distributions for those processes.Besides, it is worth noting that the technique used to estimate the warranty costs, and the subsequent computation of the simultaneous condence bands is non-parametric, and therefore no distributional assumption about the failure/repair processes is required by the proposed method.
Revista Colombiana de Estadística 41 (2018) 130 The proposed computation of the simultaneous condence bands is valid in a wide range of models that satisfy with the general conditions identied in Section 2.
The proposed computation of the simultaneous condence bands for the expected discounted warranty cost of coherent systems under minimal repair is easy to implement in current statistical package, since it only involves random sampling with replacement.Also, in the simulation scenarios studied in Section 4, reasonable actual coverage probabilities were obtained, particularly when it there was a number of systems under warranty greater than 100, a ne partition of the warranty term and a large number of resamples.
Revista Colombiana de Estadística 41 (2018) 130 This paper has focused on the minimal repair model and on the warranty policy in which repair cost is shared by both the customer and the manufacturer (PRW policy).In the literature there are references to other types of repair besides minimal repair, such as: perfect repair, in which the product is as good as new after repair, and imperfect repair, in which the product is better than used, but not as good as new.Both types of repair make use of well-known renewal processes (Park & Pham 2010, Su & Shen 2012).In addition, it is possible to con-

Proof of condition C-2
Consider δ = 1.We want to show that Using Proposition A.1 on each integral on the right hand side of the previous inequality, it is observed that there is κ 1 > 0 such that for every i = 1, . . ., m it holds that s On the other hand, Because H i (s) with s ∈ [0, t] are non-negative, bounded and continuous functions in [0, t], there is 0 < υ < ∞ such that H i (s) ≤ υ for every i = 1, . . ., m and since g (x) = x 3 with x ≥ 0, is an increasing function, from (B12), it follows that In addition using inequality ( Applying the expected value it gives } is a Non-Homogeneous Poisson Process (NHPP).Therefore, for each realization (j), it can be shown that N Revista Colombiana de Estadística 41 (2018) 130 where S x,3 are the Stirling numbers (Riordan 1937) and s 0 λ i (u) du < ∞ for each i = 1, . . ., m.Then, by applying Proposition A.1 it follows that there is κ 4 > 0 such that s x . (B16) Note that, in (B16) the sum on the right side of the inequality is a non-negative nite number and therefore there is κ 5 > 0 such that for each i = 1, . . ., m, it holds that

Appendix C. Results of Simulation Study
Some results from the simulation study described in Section 4 are showed now.

0≤s≤t|Q
(s)|, to which the process sup 0≤s≤t |U •n (s)| weakly converges is unknown.Using Corollary 1 This section uses simulation to assesses the properties of the proposed simultaneous condence bands.The simulation study considers dierent scenarios which depend on factors that may aect the performance of the bands established in Theorem 6. Revista Colombiana de Estadística 41 (2018) 130 Condence Bands for the Expected Discounted Warranty Cost of a Coherent System 13 4.1.Simulation Factors and Parameters

•
Minimal repair cost by component.It corresponds to the minimal repair value in the i-th component and is denoted c i .For this study, c 1 = c 2 = 3 and c 3 = c 4 = 5 were considered.Revista Colombiana de Estadística 41 (2018) 130 • Nominal coverage probability.Denoted by (1 − γ).It species the expected probability that the true mean cost function is bounded.A value of (1 − γ) = 0.95 was considered.• Number of simulations.This is denoted by N .A total of N = 10000 simulations were used.The Weibull distribution for component lifetime was chosen for its frequent use in industrial reliability.Besides, the values of Weibull distribution parameters set for each component ensure the distributions have increasing failure rates, and a record of failures with at least one event during the established warranty term.The constrains F 1 ≡ F 2 and F 3 ≡ F 4 were used to simplify the simulation scenarios.

Figure 1 :
Figure 1: Approximate B * (t) for systems under study.The solid curve is calculated for Hi(t) = ci exp (−t) and the dashed cure is calculated for Hi(t) =

Table 1 :
Simulation Factors and Their Levels.

Table 2
summarizes the xed values for the parameters considered in the sim-

Table 2 :
Fixed Values of Simulation Parameters.
Condence Bands for the Expected Discounted Warranty Cost of a Coherent System 15 c) Obtain q γ (t), the (1 − γ) quantile of the estimated approximate distri- j ) for each r, where A is a partition of size M of [0, W ]. Revista Colombiana de Estadística 41 (2018) 130 (Chien 2008, Chien 2010Policies (FRW) or combination of FRW and PRW policies Revista Colombiana de Estadística 41 (2018) 130Condence Bands for the Expected Discounted Warranty Cost of a Coherent System 21(Chien 2008, Chien 2010).Generalizing the results described here into previously described situations is considered relevant.

Table C1 :
Actual Coverage Probabilities for the 1-out-of-3 Components System with Hi (t) = ci exp (−t).

Table C3 :
Actual Coverage Probabilities for the 2-out-of-4 Components System with Hi (t) = ci exp (−t).