Reliability-based forecasts can be made from field data on complaints, failures, repairs, age-replacements (life limits), NTFs (no trouble found), WEAP (warranty expiration anticipation phenomenon), spares, warranty claims, or deaths. Some spares inventory forecasting software says… “Please enter forecast______” No kidding. 1800 years ago Roman Jurist Ulpian made actuarial pension cost forecasts for retiring Roman Legionnaires. Would you like actuarial forecasts? Their distributions? Stock recommendations?

## Forecast Warranty Reserves?

I asked Mike, the Apple Computer Service Department accountant, if he wanted reliability-based forecasts for estimating warranty reserves based on age specific field reliability of products and parts? He replied, “No thanks, we do that.” I asked, “What if your warranty reserve estimates are wrong?” He explained, “Not a problem, we take a charge or credit based on the error.” I wondered about better use for $$$ squirreled away in warranty reserves? I wondered how much of warranty reserves disappeared?

Mike, a tax accountant wrote, “The warranty reserve is a liability on the balance sheet so reduces the net worth of the company. For book purposes, it is also a deduction against income but, alas, not for tax purposes. Warranty expense is a tax deduction when an actual expense is paid such as providing a refund or replacement product, etc. For book purposes, when warranty expense is recognized, the reserve is reduced and that amount taken into income. The only enterprise that I can think of that can set up a reserve and deduct it is an insurance company. We have an *actuary* compute the expected payout for the current and future years.” [Struzzieri and Hussian]

Insurance companies expense actuarial forecasts of warranty reserves in year of sale, or, if not expensed, they avoid tax if claims are less than retained earnings. Unpaid-out warranty reserves (retained earnings) revert to stockholder equity, untaxed.

I thought, why not make actuarial warranty reserve forecasts and obtain the same tax advantages as insurance companies? I asked a lawyer if any companies had actually thought of doing this. She replied, “Oh, yes! That’s why Warren Buffet bought GEICO.” Which would you trust more: a forecast made by your own company reliability people or by the insurance company’s actuaries?

## Forecasting and Inventory Software

“The advantages of a simple model are:

It is economical of time and thought

It can be understood readily by the decision maker

If necessary, the model can be modified quickly and effectively”

[Bierman, Bonini and Housman, *Introduction to Quantitative Business Decisions,* Addison-Wesley, New York, 1973]. Do you want a simple forecast?

The simplest forecast of spare parts’ requirements is to guess. (Carol told me that’s how Apple Computer did it before 1990.) Less simple is the naïve forecast N*AFR = Installed base*Annual Failure Rate (or N/MTBF); this naïve forecast could be improved by accounting for attrition in installed base N; forecast auto parts’ demands as N(t)*AFR(t) times part’s demand rate AFR(t) for parts used in surviving vehicles of age t [www.SmartCorp.com and www.Epicor.com]. This requires stores to track parts’ sales by vehicle year, make, model, and perhaps engine. Does your local auto parts store collect customer’s vehicle information?

Actuarial forecasts are the SUM[a(s)*N(t-s); s=1,2,…,t] where a(s) is the age-specific (actuarial) failure rate estimate and N(t-s) is the number of retirees at each age t-s [Ulpianus 223 ad, Halley 1693, insurance 1880?]. Actuarial forecasts are the expected value of random future deaths. What if the forecast is for parts of repairable systems that continue to age after part replacement?

The actuarial demand forecast for repairable products or parts is SUM[d(s)N(t-s); s=1,2,…,t], where d(s) is the demand rate for N(t-s) parts’ installed in products of age t-s containing the parts (railroads (1920?), US AFLC (1960), Apple Computer in the early 1990s, www.Triad.com in the mid-1990s). Age-specific demand rate d(s) is not the same as actuarial rate a(s), because d(s) is conditional on the age of the repairable product containing the part(s), and a(s) is conditional on the age of the part.

Table 1. Alternative naive and actuarial forecasts with or without life data (WOLD). Nelson-Aalen or Kaplan-Meier are nonparametric estimators of cumulative failure rate function or reliability function. Either yields the same actuarial failure rate estimates.

Naïve: N*AFR | Naïve: Survivors*AFR | Actuarial Nelson-Aalen | Actuarial WOLD |

N = Installed base, population or sales in each and prior periods | N = Survivors after removing failures from installed base | Or Kaplan-Meier, with grouped life data; e.g., Nevada table | Without life data, using ships and returns counts |

ABC, MRO, and ERP software include inventory management methods. See https://www.softwareadvice.com/inventory-management/ or search Internet for “spare parts inventory”. MRO is useful for computing parts’ installed base by age from product installed base. “Gozinto Theory and Parts’ Installed base,” https://accendoreliability.com/?s=gozinto/.

Seasonality matters. Some software shows whatever historic info the customer wants: inventory turns, moving averages, lead times, what-ifs?, etc. (www.smartcorp.com and www.Epicor.com). Some software does statistics for you; e.g., “Probabilistic forecast” [www.lokad.com “No statistics required”]. Some do time-series extrapolations [www.AutoBox.com, SAS, JMP,…]; some extrapolate winner of “tournament” competition among alternative forecast models. Some dump CMMS data into Artificial Intelligence programs [Posadas et al., (LMI, United States Patent 8600843), Resienkiewicz, www.logility.com, www.lokad.com, https://www.palantir.com/impact/inventory-optimization/].

Hans Levenbach uses an “Alphabet” forecast, which compares forecast vs. demand over several periods using Kullback-Leibler divergence, the ratio of (logarithms of) probabilities. Levenbach compares f(s)=forecast(s)/SUM[forecast(s), s=1,2,…,t] and a(s)=demand(s)/SUM[demand(s), s=1,2,…,t]. Levenbach’s K-L divergence = SUM[f(s)ln(f(s)/a(s))]. K-L Divergence is intended to compare information in distributions. “Alphabet” is an empirical distribution, not the sample distribution of estimates.

People recognize the randomness of demands. Some use the Newsboy Problem solution (https://en.wikipedia.org/wiki/Newsvendor_model/) for a single period. My ex-wife’s father used the Wilson square-root formula for economic order quantity for Pep Boys circa 1960. (s,S) Inventory models provide both the reorder point “s” and the order-up-to quantity “S” [Iglehart, Vienott, (https://patents.justia.com/inventor/thomas-reed-willemain/)].

What about simultaneous parts’ inventory? Demands for different parts’ could be dependent on each other in addition to dependence on age. “Customers who bought part X also bought part Y.” Kits are collections of parts required for major repairs or overhauls such as engine or automatic transmission overhauls. What parts should be included in the kit?

At the UC Berkeley Statistics Department, Professor Bin Yu’s AI team made US county COVID-19 case forecasts that accounted for correlations among neighboring counties, to help allocate health-care resources. The same should be true of forecasts for parts. Before tackling that, let me show some actuarial single part forecasting and inventory recommendations.

## Compare Alternative Forecasts’ SSE

Tables 3 and 4 compare sums of squared errors (SSE) of failure forecasts of cohorts starting in each of 10 successive periods. (This in known in biostatistics as a staggered start model.) I.e. first cohort is observed for 10 periods, second cohort is observed for the last 9 periods, etc.: computed MTBF = 1/AFR and nonparametric failure rate estimates form simulated failures, and computed sums of squared errors in failure counts, SSE=SUM[(observed−expected)2; 1,2,…,20], where observed is simulated failure counts and expected is the forecast for four alternative forecast methods. The actuarial Nelson-Aalen or Kaplan-Meier forecast uses the grouped life data (Nevada table (yellow) in rows 1-4 of table 2), and the actuarial WOLD forecast uses the 20 unit cohorts shipped per period and the sums of failure counts grouped in each of 10 successive periods, without any record of which cohort the failures came from (bottom row of table 2).

Table 2. Example ships and failure data for a nonrepairable product or part. Yellow is “Nevada” table of grouped failure counts. Bottom row is total returns in each period.

Period | Ships | Fail 1 | Fail 2 | Fail 3 | Fail 4 |

1 | 20 | 1 | 3 | 2 | 1 |

2 | 20 | 0 | 1 | 2 | |

3 | 20 | 1 | 0 | ||

4 | 20 | 1 | |||

Sums | 60 | 1 | 3 | 4 | 4 |

Table 3. Statistics of sums of forecasts’ squared errors (SSE) from 20 simulation runs: Weibull reliability exp[(t/10)^{0.5}] with mean=20. Table lists average, standard deviations, etc. and CV = standard deviation/average. WOLD means “without life data” using the last row of sums in table 2.

Naïve N*AFR | Naïve Survivors*AFR | Actuarial Nelson-Aalen | Actuarial WOLD | |

Average | 101.82 | 100.47 | 75.13 | 164.44 |

Stdev | 34.12 | 31.27 | 34.44 | 53.68 |

Min | 44.79 | 49.17 | 26.24 | 66.52 |

Max | 161.15 | 156.57 | 139.01 | 291.55 |

CV | 33.51% | 31.13% | 45.85% | 32.64% |

Table 4. Statistics of sums of forecasts’ squared errors from 20 simulated runs: Weibull exp[-(t/10)^{1.5}] with mean=9.03

Naïve N*AFR | Naïve Survivors*AFR | Actuarial Nelson-Aalen | Actuarial WOLD | |

Average | 71.31 | 71.58 | 67.25 | 256.48 |

Stdev | 47.78 | 44.71 | 52.59 | 94.08 |

Min | 19.98 | 22.63 | 18.00 | 126.50 |

Max | 236.20 | 224.76 | 244.31 | 520.11 |

CV | 67.00% | 62.46% | 78.21% | 36.68% |

Naïve vs. actuarial with grouped life data (Nelson-Aalen or Kaplan-Meier) averages less SSE and comparable standard deviation of SSE vs. without life data. Surprisingly naïve and actuarial with grouped life data have larger coefficients of variation of SSE vs. without life data, despite smaller averages and standard deviations of SSE. Without life data gives more scatter among hindcasts, due to Less information in period returns counts than in life data.

## What is the distribution of demand?

We need to forecast the distribution of demand to set order quantities, stock levels, and reorder points. The probability an order is filled upon demand is the “fill rate.” The probability of a stock-out between inventory refills is the “service level.” How to manage inventory depends on your objective(s). Minimize backorders? Their durations? Minimize holding costs subject to service level or fill rate and space constraints? Maximize profit from service operations? All of them? They all depend on conditional demand rate d(s). The demand distribution forecast should include the sample uncertainty in actuarial demand rate estimates!

Forecasts change as installed base and its ages change. Recommendations vary: [King, Willemain et al., Veinott, Iglehart,…]:

- Percentile of Poisson[forecast] (US Air Force Logistics Command, others)
- Percentile of normal distribution N[forecast, variance]. Some people use the actuarial forecast as the mean and guess the variance.
- Percentile of negative binomial distribution [patent number 5287267 by Robin Roundy]
- Upper prediction limit of normal distribution N[forecast, variance] with estimated variance
- Upper prediction limit of estimated non-normal demand distribution: e.g. for frequent zero-demand periods [Ghobbar and Friend].
- Plug into K-L divergence [www.lokad.com, Hans Levenbach]

In statistics, a tolerance limit or a prediction limit both put bounds on variation in the future. Theprediction limit bounds a single future sample period, whereas a tolerance limit bounds the entire population (equivalently, an arbitrary sequence of future sample periods). [Wikipedia] Compute the prediction limit for a demand forecast, if you’re forecasting one period ahead. Compute tolerance limits if you’re trying to forecast the entire future [Hall and Wellner, Wang, Jewell, and Tsai].

The simulation in tables 3 and 4 was for cohorts of 20 units put into service in successive periods. You could track each unit from birth to death, group ages at failures from each cohort by calendar period, or you could use ships and failures or returns counts grouped, without regard to cohort of origin, using data required by generally accepted accounting principles.

It is unusual to have the same ships every period as in table 2. Ships and returns counts can be modeled by a nonstationary Poisson-Exponential M/G/infinity self-service system, where M stands for Poisson distribution of periodic arrivals or ships cohorts (or M(t) for nonstationary Poisson), G stands for the probability distribution of service time or life time, and infinity means that every arrival serves for a random service time with reliability function 1-G(t). It is difficult to disprove ships cohorts have nonstationary Poisson distributions from a single realization of staggered-start ships cohorts.

The nonparametric maximum likelihood and least squares estimators estimate actuarial failure rates a(s) [George and Agrawal 1973, Mirasol, Newell] and age-specific demand rates d(s) [George]. The joint distribution of actuarial rate estimates a(s) or d(s), s=1,2,…,t is asymptotically multivariate normal according to the martingale central limit theorem [George 2021]. The mean vector and variance-covariance matrix can be computed or obtained by simulation, resampling (bootstrap), or from the Cramer-Rao limits (I use Mathematica to compute the inverse of the Fisher information matrix.) The actuarial forecast is a sum SUM[d(s)n(t-s)] of multivariate normally distributed demand rates d(s) times the installed base n(t-s), so it is reasonable to expect the forecast to also be asymptotically normally distributed, if the installed base n(t-s) is known and not random. (This is not necessarily a contradiction to the Poisson output of and M/G/Infinity self-service system, because the normal distribution is an approximation to the Poisson distribution under some conditions.)

The demand forecast distribution depends on joint distribution of the d(s), s=1,2,…,t estimates. Assuming the installed base n(t-s) is known, the variance of the demand forecast is (Notice the covariance?)

Var[SUM[d(s)n(t-s)]] = SUM[Var[d(s)]n^{2}(t-s)]+SUM[Covar[d(s),d(u)]n(t-s)n(t-u), s, u = 1,2,…,t].

Table 5 shows results of 2000 simulations of a two-period, two staggered-start cohorts of Poisson ships at rate 20 per period with probabilities of failure 0.1 and 0.05 at ages one and two periods. The left block of four columns contains the statistics from the Kaplan-Meier or Nelson-Aalen estimator from life data, and the right block of four columns contains the statistics from the nonparametric estimator from ships and returns counts, without life data. (Both are maximum likelihood estimates.) The actuarial rate estimates from the Kaplan-Meier estimator are more accurate and precise. The actuarial rate estimates from the Kaplan-Meier estimator have zero correlation theoretically. The actuarial rate estimates without life data are correlated so this must be accounted for in the variance of forecasts.

The martingale central limit theorem says actuarial rate estimates from ships and returns counts has asymptotic multivariate normal distribution, *with correlations.* That implies the demand distribution also has asymptotic normal distribution. The skew and kurtosis of a1 and a2 actuarial rate estimates are constrained by the facts that actuarial rates are nonnegative and less than or equal to 1.0. Is the accuracy and precision worth the cost of tracking products and parts for life data?

Table 5. Compare alternative nonparametric estimates of actuarial demand rates a1 and a2 from grouped life data (columns 2-5) vs. ships and returns counts (columns 6-9).

mean a1 | mean a2 | var-cov | var-cov | mean a1 | mean a2 | var-cov | var-cov | |

Mean | 0.100601 | 0.050225 | 0.002489 | -4.3003E-05 | 0.089014 | 0.075998 | 0.003501 | -0.00206904 |

Stdev | 0.049886 | 0.052545 | -4.3E-05 | 0.002760965 | 0.059173 | 0.088822 | -0.00207 | 0.00788927 |

Skew | 0.469606 | 0.919846 | correl | Correl | 0.445344 | 1.335437 | correl | Correl |

Kurt | 0.246481 | 0.444015 | 1 | -1.64% | -0.09186 | 1.873055 | 1 | -39.37% |

CV | 0.495878 | 1.046193 | -1.64% | 1 | 0.66476 | 1.168741 | -39.37% | 1 |

## More Work

Forecasting usually means forecasting the future, at least one period beyond the known ships and returns periods. That requires estimating d(t+1) and the next cohort of ships. I usually extrapolate a(t+1) or d(t+1) by regression, which gives standard errors of the extrapolations. Same for the next cohort’s ships.

Multiple product or parts’ optimal opportunistic replacement call for estimates of the joint distributions of parts’ actuarial demand rates. What else should we replace as long as we’re replacing one part? Kits? Some people compute distributions using regression extrapolation of actuarial rates. Some people use bootstrap, jackknife, resampling, or other simulation. Multivariate: demands are correlated; e.g., copulas! Multi-Echelon inventory complicates forecasting; e.g., US AFLC: engines go from bases to Tinker or Kelly AFB for overhaul and back to bases, depending on flying-hour program, data, and other circumstances. Send data and challenges; then maybe I’ll have more to tell you.

## References

George, L. L. and A. Agrawal, “Estimation of a hidden service distribution of an M/G/∞ system,” *Naval Research Logistics,* 20: 549–555. doi: 10.1002/nav.3800200314 “M/G/Infinity Service Distribution,” https://sites.google.com/site/fieldreliability/would-you-believe-you-dont-need-life-data/mginfinity-service-distribution/ 1973

George, L. L., “Renewal Process Estimation, Without Life Data,” https://accendoreliability.com/?s=renewal/ 2021

George, L. L., “Actuarial Forecasts, Least Squares Reliability, and Martingales,” https://accendoreliability.com/actuarial-forecasts-least-squares-reliability-martingales/#more-421021/, 2021

Adel A. Ghobbar and Chris H. Friend, “Evaluation of forecasting methods for intermittent parts demand in the field of aviation: a predictive model,” City University London, School of Engineering, *Computers & Operations Research* 30, 2097 – 2114, 2003

Hall, W. J. and Wellner, J. A. “Confidence bands for a survival curve from censored data” *Biometrika* 67, 133-43, 1980

Iglehart, D. L. “Optimality of (s, S) Policies in the Infinite Horizon Dynamic Inventory Problem” *Mgmt. Sci.* 9, 259-267, 1963

Peter L. King, “Understanding the Safety Stock and Mastering its Equations,” *APICS Magazine,* July-August 2011

Hans Levenbach, “Why Demand Planners Need to Improve Data Quality Before Using Forecasting Models for Planning and Budgeting,” https://www.linkedin.com/pulse/how-select-best-performing-forecasting-methods-supply-hans/, Dec. 2020

Noel M. Mirasol, Letter to the Editor—”The Output of an M/G/Infinity Queuing System is Poisson” *Operations Research* 11(2):282-284. https://doi.org/10.1287/opre.11.2.282, 1963

G. F. Newell, “The M/G/Infinity Queue,**” **SIAM *Journal on Applied Mathematics,* Vol. 14, No. 1, pp. 86-88, Jan., 1966

Sergio Posadas, Carl M. Kruger, Catherine M. Beazley, Russell S. Salley, John A. Stephenson, Esther C. Thron, Justin D. Ward, “Forecasting Parts Demand Using Service Data and Machine Learning,” Volume 1. Logistics Management Institute, Jan. 2020

Maria Rosienkiewicz. “Artificial Intelligence Methods in Spare Parts Demand Forecasting.” Wrocław University of Technology, Poland, https://www.researchgate.net/publication/298203927, Jan. 2013

Paul J. Struzzieri, FCAS Paul R. Hussian, FCAS, “Using Best Practices to Determine a Best Reserve Estimate,” https://www.casact.org/sites/default/files/database/forum_98fforum_struhuss.pdf/

Veinott, A. F., Jr. “The Status of Mathematical Inventory Theory” *Mgmt. Sci.* 12, 745-777, 1966b

Mei-Cheng Wang, Nicholas P. Jewell, and Wei-Yann Tsai, “Asymptotic Properties of the Product-Limit Estimate Under Random Truncation,” *The Annals of Statistics,* Vol. 14, No. 4, pp. 1597-1605, 1986

Thomas R. Willemain, Charles N. Smart, Henry F. Schwarz, “A new approach to forecasting intermittent demand for service parts inventories”, *International Journal of Forecasting* 20, 375 – 387, 2004

Thomas Willemain, “Top 3 Most Common Inventory Control Policies,”

https://smartcorp.com/inventory-control/inventory-control-policies-software/, 2019

Larry George says

Wonder where the Yellow went? Table 2 “Yellow” refers to the grouped failure counts in periods 1-4, the triangular set of numbers looking like Nevada on its side.

I forgot the minus sign in Weibull exp[-(t/10)^0.5]. It’s OK in my spreadsheet.

Table 5 was too wide so the important info didn’t show. The visible left part is statistics of actuarial rate estimates from grouped life data. The invisible part is statistics without life data, just ships and returns. Notice that the parameters without life data are a little worse. Better than nothing! Contact me if you want the original *.doc or *.pdf.

avg a1 avg a2 var-cov avg a1 avg a2 var-cov

Mean 0.101 0.050 0.002 -4.3E-05 0.089 0.076 0.004 -0.002

Stdev 0.050 0.053 -4.3E-05 0.003 0.059 0.089 -0.002 0.008

Skew 0.470 0.920 Correl 0.445 1.335 Correl

Kurt 0.246 0.444 1 -1.64% -0.092 1.873 1 -39.4%

CV 0.496 1.046 -1.64% 1 0.665 1.169 -39.4% 1