Forecasting the returns of reusable containers

JOURNAL OF OPERATIONS MANAGEMENT Vol. 8, No. I, January 1989

Forecasting the Returns of Reusable Containers PETER KELLE* EDWARD A. SILVER**

EXECUTIVE SUMMARY A number of organizations sell products in containers that can be reused. Illustrative cases include beverages (in kegs, plastic cartons, etc.) and liquid gases (in cylinders). The time from issue to return of an individual container is usually not known with certainty and there is a chance that the container is never returned (because of loss or irreparable damage). In such a situation, even if the demand pattern is known and level with time, it is still necessary to acquire new containers from time to time. Such an acquisition must be initiated when the inventory level of containers drops too low in relation to the anticipated net demand (demand minus returns) during the replenishment lead time. Thus, it is important to forecast the net demand as well as to obtain an estimate of the accuracy of the forecast. In this paper, four different forecasting procedures, based upon different amounts of information, are developed. The information possibilities include recent issues period by period and recent returns (either individually identified or on an aggregate basis). The methods are compared on a wide range of simulated data, including some cases based on empirical data obtained from industry. Not surprisingly, the use of additional information improves the performance. However, most of the benefit associated with using the costly approach of identifying and tracking the issues and returns of individual containers is achieved by the more practical method of recording only aggregate issues and aggregate returns period by period. INTRODUCTION A n u m b e r o f organizations sell products in containers that can be reused. An adequate supply o f e m p t y containers must be on hand to satisfy the demand. Part o f this supply is a result o f the returns o f previously issued containers. The time from issue to return o f an individual container is a r a n d o m variable with a distribution that includes a finite probability of never being returned (lost, d a m a g e d , etc.) In order to properly establish the reorder point for purchasing new containers, it is necessary to forecast the net d e m a n d (demand minus returns) during the purchasing lead time and the variability o f this net demand. A complicating factor is that during the lead time the s a m e container may be issued, returned, reissued, etc. In particular, the returns late in the lead time can d e p e n d upon the issues (demands) earlier in the lead time and such issues are t h e m s e l v e s r a n d o m variables w h o s e values are typically not known at the start o f the lead time. The problem studied in this paper was motivated by interactions (visits and consultations) with organizations selling products in returnable containers. These include i) liquid gases (in

*Louisiana State University, Baton Rouge, Louisiana **University of Calgary, Calgary, Alberta, Canada

Journal of Operations Management

17

cylinders), ii) beer (in kegs), and iii) non-alcoholic beverages (in returnable bottles and/or plastic cartons). Although detailed data on demands and returns were not readily available from these sources, approximate parameter values were obtained for use in our simulation experiments to be discussed later in the paper. Some authors have dealt with the forecasting of returnable containers, mainly concentrating on the estimation of the average number of times that an individual, container is reused before it is lost or irreparably damaged (see Flies (1977), Fujita (1977), and Weymes (1978)). Recently, a Box-Jenkins statistical approach was applied to the estimation of the return distribution by Goh and Varaprasad (1986). None of the above papers has focused on the lead time return and net demand forecast taking account of the correlated nature of the random variables and the effects of random demand during the lead time itself. The estimation of a measure of forecast errors (e.g., variance or mean absolute deviation) was also not considered in the earlier studies. In this paper different forecasting methods, dependent upon the available data, are developed to estimate the returns and net demand during the lead time. A measure of forecast errors and the appropriate reorder point are also estimated for each of the forecasting methods. The next section of the paper describes four different forecasting methods, each one based on different available information. Method 1 is the simplest case. It utilizes only • the expected value and the variance of the demand during the lead time, and • the probability of each container eventually being returned. Method 2 uses more detailed information, namely • the actual issues during each previous period and • the probability of return in 1, 2 ..... n periods for any given container. Method 3 uses the same issue and return probability data as Method 2, as well as • the amount returned up to the present from each previous issue. This additional information permits us to obtain the appropriate conditional probabilities of return quantities during the lead time for the remaining outstanding portion of each previous issue. As an extreme illustration suppose that an issue of 20 containers was made in period 2, we are currently at the end of period 4, and all 20 containers have already been returned. Then, clearly the conditional probability of 0 returns from this particular issue during the lead time must now be 1. These conditional probabilities permit us to calculate the best possible forecast based on past data. Unfortunately, to acquire the additional information, the containers have to be individually identified and records must be kept as to when they were issued and returned. To avoid this substantial cost and still partially utilize information on previous returns, we developed the following fourth method. Method 4 requires, besides the issue and return probability data of Method 2, only • the total amounts returned in each of the recent periods without identification of when the associated containers were issued (i.e., only aggregate return data). The next section of the paper evaluates the four different forecasting methods. The comparison is made by means of a wide range of simulated data including some empirical data from industry. Finally, some conclusions are drawn with respect to the use of the different forecasting methods.

18

FORECASTING METHODS Method 1: Forecast Based on Average Behavior

In this simple forecasting method we assume, as is sometimes done in practice, that each demand has a fixed probability P of an accompanying return of a reusable container. The only pieces of information that are utilized are i) the probability P that a container is ever returned, and ii) the expected value E(dE) and variance Vat(dE) of the random lead time demand dE. Under the above assumptions, for a given value of dE, the lead time return is a binomial random variable since each of the dE demand units has the same probability P of an accompanying return. A complicating factor is that dL, the first parameter of the binomial distribution, is itself a random variable. The distribution of the resulting so-called mixed binomial random variable can be expressed in a closed algebraic form only for Poisson distributed dE, which can be a poor approximation if the empirical mean and variance of dE are not close to each other. In Equations (A. 1) and (A.2) of Appendix A, we have expressed the exact mean and variance of the mixed binomial in simple closed algebraic forms for arbitrary random dE with known mean E(dE) and variance Var(dL). Hence, the expected lead time returns under Method I, denoted by ER Cll, can be expressed as ER CI) = P E(dL).

(1)

This quantity is then used as a forecast of the lead time returns. The variance VR m of lead time returns has the form VR¢I) = pZ Vat(dr.) + P(I-P) E(dE).

(2)

In deciding on the need to acquire new containers the net demand (demand minus returns) during the replenishment lead time has to be forecasted. The expected lead time net demand (the forecast) using Method 1, and denoted by ED (~, is simply EDtt) = E(dL) - ER li) = (l-P) E(dE).

(3)

To summanze, the expected lead time net demand is the expected demand, E(dL), minus the expected returns, the latter being the product of P and E(dL). In the calculation of the variance of lead time net demand, VD, we have to account for the correlation between the random dL and the random lead time return. Using the algebraic expression (A.3) derived in Appendix A, we have VD Cl~ = (I-P) 2 Var(dE) + P(I-P)E(dL).

(4)

The reorder point s is expressed in the common way used in the inventory control literature (see e.g., Silver and Peterson (1985, Chapter 7) for details), as s = ED + k* ~

(5)

where k* is the appropriate safety factor (based on service considerations or on minimizing expected total relevant costs). The forecasting method is illustrated on a numerical example in Table 1. The table also shows illustrations of the other methods to be described in the next three subsections.

Journal of Operations Management

19

TABLE 1 NUMERICAL EXAMPLE ILLUSTRATING THE FOUR DIFFERENT FORECASTING METHODS Return distribution: Pl = 0.4, p,.. = 0.2, P3 0.2; Demand parameters: p, = 30, o = 6; lead time L = 2; actual period: t = 5; k* = 2.32. =

t

OBSERVED VALUES OF ISSUES AND RETURNS: ?ERIOD

ISSUE

RETURN 1 PERIOD LATER

RETURN 2 PERIODS LATER

RETURN 3 PERIODS LATER

1

34

12

8

8

2

22

12

3

4

3

35

10

5

4

35

12

5

31

FORECASTS: METHOD

I Based on aggregate behavnor

2 Based on recenl issues

INFORMATION USED

RETURN FORECAST

NET DEMAND

P = 0.8

Ill ERClI=0 8x60=48

131 E D ' " = 12

E(dL)=60 Var(dL)= 72

[21 VRII)=0 82x72+0 8x0.2x60~55 68

141VD I l l = 0 2-'x72 + 0 8x0 2x60 = 12.48

Pl = 0 4 p 2 = 0 2 p a = 0 2

161 E(WO= 35x0.2= 7

[131 EDt2*=60-51 6 = 8 4

u~= 35 u4 = 35 u~= 31

17] Var(WO= 35x0 2x0 8 = 5 6

1141 VDC21=6" + 5 6 + 8 4 + 7 4 4 +

REORDER LEVEL

[51 s " l ~ 2 0

sq:)= 30

only E(W4) = 14 E(W~)= 18.6

(6x0.8) 2 + 30x0 2x0 8

Vat(W4)= 8 4 Var(W~) = 7 44

= 87 28

[9] E(W6)= 30x0 4 = 12 1101 Vat, W,) = (6x0 212 + 30x0 2x0.8 = 6 24 1111ERq2)=51 6

[121VR~21= 27 68 1181 ED 131= 60-55.93 = 4.07

3. Based on issues and individual udentlfied return.*,

Pl = 0 4 p 2 = 0 2p3 = 0 2

[ 161 ER °1 = 55 93

117] VR O)= 23 79

4

Pu = 0 4 p 2 = 0 2 p a = 0 2

1221 E C = 3 3

1241VC=3 9

126] EDI4J=5 I

ua=35 ua= 35 u~= 31

[231 ER q41= 54 9

1251 VR c41= 23 78

(27) VDI41=8I 9

Based on issues and aggregate relurlls

[151 E(W~[ V3)=(35-15)x0.~0 4 = 10

ua= 3 5 ~ = 35 us= 31

Var(W~ I V31 = 20x0 5x0.5 = 5

V~= 15 Va = 12

E(Wa [ V.O= 15 33 Vat(W41V41 = 5 11

s~31= 25

1191VDq3~=81 8

s~4)= 26

y4 = 21 ys= 21

*Note: A bracketed ([ ]) number refers to the number of the equation used. M e t h o d 2: Forecast B a s e d O n l y on R e c e n t Issues In practical situations usually more information is available than was used in Method 1. We n o w develop a more exact forecasting m e t h o d based on the f o l l o w i n g data: i) the actual issues u~ during each previous period i <__ t where t denotes the last period that w e have observed (i.e., w e are currently at the end of period t), and

20

ii) the probability distribution of the number of periods from issue to return of any single container, pj denotes the probability of return after exactly j periods. The largest j for which pj is significant (non-zero) is denoted b y n (n takes on values in the range of 2 to 20 time periods in practice). The sum x pj = P is the total probability of a J=l

a return, used in Method 1 and p~ = 1-P is the probability of a container never being returned. The return distribution (Pl,P2 ..... p,,p~) is an arbitrary discrete distribution that can be specified from observations of returns of individually identified containers. If such data are not available, statistical methods, such as regression analysis or the identification and estimation phases of the Box-Jenkins procedure, the latter as described by Goh and Varaprasad (1986), can be utilized to estimate the return distribution. It should be noted that any time needed for inspection, cleaning and repair of returned containers can be included in the time from issue to availability for use rather than just from issue to return. In our analysis we assume that the probability of the return of a particular unit is independent of the retum of other units. Hence, the returns from the issues (the ui's) are a sequence of independent multinomial trials with the success probabilities pj (j = 1..... n). The returns in different periods from a particular issue ui are definitely correlated, but the returns from different issues are independent. (The mathematical details of the multinomial model are described in Appendix B.) Suppose that we are currently at the end of period t and the lead time is of duration L periods. Then the total return Yt.L during the lead time periods t + 1 ..... t + L is the sum of the independent lead time returns Wi from the issues ui of the different periods i < t + L-1. Each W, is a binomial random variable with a number of trials u, and success probability Ri of each particular unit returning in one of the lead time periods t + 1 ..... t + L . Ri = Ri(Pj'S) is a function of i and the return distribution. It is a partial sum of the pj's as expressed in (B. 1) of Appendix B. We now treat previous issues (i <__ t) and issues during the lead time (i > t) separately. The previous issues ui (i <__ t) are observed, known values, thus the total lead time return from the previous issues is the sum of independent binomials with known parameters. Its expected value is ER~ 2) =

t t ~ E(Wi) = ~. i= 1 i=im

ui Ri(pj's)

(6)

and variance t

VR[ 2) =

~.

t

Var(W0 =

i= 1

with

Y-, i=im f

u~ R~(pfs)

/

im = m a x ~ l , t - n + l j "

[1 - R~(pj's)]

(7)

(8)

since for i < im E(Wi) = 0 and Var(Wi) = 0. For known lead time issues the expected value and variance of the returns during the lead time from issues made during the lead time can be expressed similarly, taking the summation over i from t + 1 to t + L-1. However, the lead time issues u~ (t < i < t + L-1) are usually not known in advance. Thus the lead time return Wi from each random issue u~ is a mixed binomial

Journal of Operations Management

21

for i > t. Each random issue can be approximated by the random demand of that particular period with sufficient accuracy if the average shortage level is negligible relative to the amount of issues. (This can be ensured, e.g., by the use of a high service level in the inventory control and purchasing policy.) Under such conditions E(u,) and Var(u,) can be replaced by the expected value Ix, and variance cri2of the demand in period i. Applying the results (A. 1) and (A.2) of Appendix A for the total lead time return of the lead time issues we get the expected value ERE =

t+L-1 Y. i=t+l

E(W0 =

t+L-1 E i=t+l

Ix, R,(pj's)

(9)

{[cr, R,(pj's)] 2 + Ix, R , ( p j ' s ) [ 1 - R , ( p j ' s ) ] }

(10)

and the variance VRL =

t+L-I t+L-1 Var(W,) = E .~ i=t+l i=t+l

where again the R,(p/s) are partial sums of the pj's defined in (B. 1) of Appendix B. Note from (7) and (10) that the variance of each retum W, from a random issue u, has the additional term [(r, Ri(pj's)] 2 compared to the case of known issue. The total lead time retum forecast by Method 2, denoted by ER (2), is the sum of expected retums of previous issues (6) and of lead time issues (9) ER C2) = ER(~ ) + ERL.

(11)

The variance of lead time returns, VR (2), is the sum of the variances (7) and (10) VR (2) = VR(t ) + VRL,

(12)

since the independence property holds for their components. The forecast of the lead time net demand can be expressed as the expected lead time minus expected lead time return ED (2) =

t+L t+L [~ E ix,- ER (2) = Y-, ~,i=t+l i--t+l i

u,R,(pj's) + 1

t+L-I ] E Ix,R,(pj's) . i=t+l

(13)

The variance of the lead time net demand is influenced by the correlation of random lead time issues and their retums. Here the result (A.3) of Appendix A can also be applied for each of the lead time issues, except ut + L which has no return in the lead time. Thus, VD (2) = O']+L + VR(t 2) +

t+L-1 E { o ' ] [ 1 - R i ( p j ' s ) ] 2 + ~ , R i ( p j ' s ) [ 1 - R i ( p j ' s ) ] } (14) i=t+l

with VR(~ ) expressed by (7) and Ri(pj's) expressed by (B. 1). Method 3: Forecast Based on Recent Issues and on Individually Identified Recent Returns

This method assumes that the containers are individually identified and records are kept of when they are issued and returned. Thus, we have the following information (the first two portions are the same as in Method 2): i) the actual issues u, during each previous period i < t, ii) the probability pj of a return after precisely j periods (j = 1,2 ..... n)

22

iii) the amount V, returned up to and including the present period t from each previous issue ui (i < t). The returns in different periods from a particular issue u, are highly correlated. Thus, using information about the realized returns of a specific issue up to a particular moment in time should permit a more accurate forecast of the subsequent returns from the same issue. In fact, this is the most accurate forecast possible using only past data. However, observation of the realized returns of particular issues requires coding of containers by issue date, recording of issue date at time of return and new coding at each issue. This type of individual identification of the returns is very rare in practice. We use this Method 3, developed for the maximal information case, primarily as a best-case benchmark for the evaluation of the other forecasting methods that require less data. We now describe the method of updating the conditional return probabilities. Recall that W, was defined as the return during the lead time (periods t + 1,t + 2 ..... t + L) from the issues ui in period i. Using the multinomial return model (see Appendix B), the conditional distribution of Wi given the observation of V, is binomial with the parameter values ui - Vi

(number of trials)

and Q,(pj's) = R,(pj's) t-i

(success probability)

(15)

1 -~pj j=l for i < t. Thus improved forecasting is possible for the lead time returns from issues u, (i < t). For the issues of periods i > t the same forecast can be applied as given by Method 2. The total lead time return forecasted by Method 3, ER ~3~, is thus t-1 ~ (u,- V,)Q,(pl's) + ut Rt(pj's) + ERL (16) i=im with im defined by (8), Rt(pj's ) defined by (B. 1) and the expected lead time retum from the lead time issues, ERL, expressed in (9). Similarly the variance of the lead time returns of issues u, for i < t can be more exactly estimated using the conditional variance of each W, given its corresponding observed Vi. The variance of the total lead time retum. VR TM, is the sum of the variances of independent binomials ER ~3) =

t-I

VR'3'=

Y_.,

r

+ u,R,(p;s)tI-R,
~

+ VRL

(17)

i =im with the same notation as (16) and with the variance of the lead time return from lead time issues, VRL, expressed by (10). The forecast of lead time net demand is similar to (13) with ER 12> replaced by ER°k Thus ED ~3) =

t+L [ t-l= t+L-I ] 1~ Ix,.~, (ui-V,)Qi(pj's) + utRt(pj's) + Y-, Ix,R,(p/s) i=t+l i im i=t+l

(18)

The variance of lead time net demand is also similar to (14) with VR~t ) replaced by VR {3) = VR 13) - VRL, since the random lead time issues and their returns are not influ-

Journal of Operations Management

23

enced by the earlier observations. The result is VD C3) = O-2+L +

t-1 E i=im

(ui-Vi) Q i ( p S s ) [ ( 1 - Q , ( p S s ) ]

+ u, R , ( p j ' s ) [ 1 - R , ( p j ' s ) ]

t+L-I

{(I .2 [l - RI(pj's)] 2 "-t'- ~,.£1Ri(pj's ) [1 - Ri(pj's)] }

+

(19)

i=t+l Method 4: Forecasting Based on Recent Issues and on Observations of Recent Aggregate Returns

In most practical situations it is very expensive or even impossible to identify and track individual containers. In this case in period k (k <_ t) only the aggregate return Yk can be observed, that is the sum of returns, that occur in period k, from all the different issues ui (i < k). Thus, Method 4 utilizes the following information (again the first two portions are as in Method 2): i) the actual issues ui during each previous period i < t, ii) the probability pj of a return of a container after exactly j periods (j = 1,2 ..... n), iii) the total amount Yk returned in each previous period (k < t). The aggregate returns, Yk (k <__ t) are correlated with the lead time returns, W f s , from the previous issues u,'s (i < t), but independent of the returns from the issues u,'s, i >t. Thus, improved forecasting is possible for the total lead time returns, i-z i w,, from previous issues by taking the conditional expected value of i= t~ i W, given the yk's (k < t). Recall that the return distribution pj has Pn as the last positive term. Thus, E(W,) = 0 and Var(W,) = 0 for i< i m = max{l,t-n+ 1}. Furthermore, the returns Yk for k < t-n+ 1 are independent of the returns Wi, i > ira. Thus, in the following random vector, W,, Yt,Yt-I..... Yt-n+2]

(20)

\i=im each component is the sum of independent binomials with different parameters. However, the components of (20) are definitely correlated. The correlation matrix can be expressed exactly on the basis of the multinomiai return model (see Appendix B). The conditional expectation EC = E

~ Wi[yt,Yt-t ..... Yt-n+2 i = im

)

(21)

cannot be expressed in an exact analytic form. However, it is reasonable to use a multidimensional normal vector (with the same variance-covariance matrix) as an approximation of the random vector (20). The normal vector should be a good approximation if the issues are reasonably large and the return distribution has several significant positive terms (n = 4 terms are enough in practice). The joint normal approximation allows us to express the conditional expectation (21) in an analytic form easily evaluated with a simple computer routine:

\i=im

24

"

with the following notation: Y = (Yt,Yt-I..... Yt-n+2) is the vector of recent aggregate return observations; T is the variance-covariance matrix of the random vector y, expressed exactly in (B.2) of Appendix B on the basis of the multinomial return model; T-J is the inverse of the matrix T; c_ is the vector, expressed by (B.3), of the covariances cI = coy

(It, \i=

Y. W,,yt_j+l im

and E ( t~1 W , ) \ i = i,ll

)

, j = l ..... n-l;

is the unconditional expectation, given by Method 2 in equation (6).

The forecast of the total lead time return by Method 4, denoted by ER C4~, is then ER (4) = EC + u, Rt(pj's) + ERL = ER t2) + c T " ( y - E ( y ) )

(23)

which shows that the forecast Method 2, based on recent issues only, is here corrected by a term depending on the difference between observed and expected values of the recent aggregate returns, an intuitively appealing result. The estimation of the variance of lead time returns of previous issues can similarly be improved using the conditional variance: (t-1 VC = Var

I ) • Wi Yt,Yt-i..... Yt-n+l \ i =im

= Var

( t - 1 ) E W, i =im /

where c ' is the transpose of vector c. Then the variance of total lead time returns, using Method 4, is VR (4) = VR (2) - c T -I c ' .

- cT -I c '

(24)

(25)

This is seen to be the estimate given by Method 2 adjusted by a term depending upon only the variance-covariance matrix of the returns. The forecast of lead time net demand has the same expression as (13) with ER (2) replaced by ER (4). Thus, ED (4) =

t+L ,S_, ixi - ER c4) = ED C21 - c__T -I (y - E(y)) i=t+l

(26)

The variance of lead time net demand is given by expression (14) of VD ~2) with VRt~) replaced b y V R ~ ) - c T -1 c I VD (4) = VD (z) - c T "1 c l

(27)

EVALUATION OF T H E F O R E C A S T I N G M E T H O D S

Evaluation Measures We used the maximal information situation (Method 3) as a basis for comparison. The performance of any one of the three other forecasting methods, relative to the base case, depends upon the actual issues and returns, both of which are random variables. Hence, the

Journal of Operations Management

25

relative performance is a random variable, and thus we need statistical measures of the relative performance. The mean deviation is not a proper measure because positive and negative errors can cancel each other's effects. We have used the following statistical performance measures to evaluate the forecasting methods: i) mean positive and mean negative deviation (to test bias), ii) mean absolute deviation or standard deviation (both common measures of variability) and iii) histogram of deviations (for detailed analysis). The Evaluation Methods

Analytic expressions, based on the probabilistic behavior of the relative performance (error), would be the most effective evaluation tool. However, the relative error (the estimated value minus true value all divided by the true value) is a complicated expression since it is the ratio of sums of correlated binomials. Neither the distribution nor even the standard deviation of the relative error could be expressed in appropriate analytic form. Approximations of the standard deviation of the relative error can be achieved by replacing the sums of binomials by normal variables. Unfortunately, the standard deviation of the ratio of correlated normal variables cannot be expressed in a simple form. In the literature different approximations are described (e.g., Hinkley (1969); Marsaglia (1965); and Hayya, Armstrong, and Gressis (1975)) but the simpler approximations are very inaccurate in our case because of the high correlation between the numerator and the denominator as well as the high coefficient of variation of the denominator. Approximations with sufficient accuracy are functions of the first four moments and are too complicated to use. Simulated random historical data were used to estimate values of the statistical measures of relative performance because the amount of data obtained from industry were not sufficient to derive statistical performance measures for a broad range of parameter values. The underlying simple random processes were generated according to specific distributions: truncated normal distributions for the demands and multinomial distributions for the returns through time from each particular issue. A comment is in order with respect to the assumed demand distribution for containers. The truncated normal is a reasonable approximation of the period by period demand for the product (issued in the container). We are effectively assuming that container requirements match the demand requirements period by period. In other words, we are ignoring the complicating effects of production (container filling) scheduling logic, including capacity considerations. In the following we describe the parameters of the simulation program and their values for which tests were executed. a) The return distribution is given by the values of pl ..... p,. The forecasting methods were tested for the following characteristics of the return distribution: n, the number of positive probabilities in the return distribution--values used were 2,3,4,5,9 and 15; n

P =

~ p3, the total probability of return, had the following j=l values: P = 0.4, 0.6, 0.75, 0.8, 0.85, 0.9, 0.96;

The shape of the return distribution may also influence the performance of the forecasting

26

methods. For that reason uniform, monotonic increasing, monotonic decreasing, unimodal, bimodal and U shaped distributions were tested. b) The lead time influence on the forecasts was examined for the values L = 1, 2, 3, 4, 5, 6, 8, 10, 20, 30. c) The demand was generated according to a normal distribution where negative values were replaced by 0. The mean values were ix = 20,30,50 and the coefficients of variation o'/1~ = 0,0.05,0.1,0.2,0.25,0.4. d) In setting the purchasing reorder point as the lead time net demand plus the safe~factor k* times the standard deviation of lead time net demand, we used the k* values shown in Table 2. TABLE 2 k* VALUES AND ASSOCIATED SERVICE LEVELS Safety factor k* Associatedchanceof a stockout during L

1.28

1.65

2.32

0.1

0.05

0.01

Each simulation run involved using parameters from a) and c) to generate time series of issues and corresponding time series of returns from each of the issues. Using this data forecasts were made and reorder points were computed by means of the four different forecasting procedures, utilizing parameters values from b) and d). To avoid problems of forecasting with inadequate history, forecasts were not done for t < 2n (the aggregate method requires 2n historical periods whereas methods 2 and 3 require only n periods). Evaluation Measures for the Simulation Runs For each combination of the parameter sets a) to d) the simulation procedure calculated and compared the forecasted values of the different forecasting methods for a large number of independent realizations of issues and returns which were generated according to the distributions specified by the input parameters. A realization here refers to a time series of demands and issues sufficiently long to permit a single lead time forecast by each of the four methods. As is common practice in the use of simulation, some trial runs were first performed to estimate the number of realizations necessary to ensure stability of the various empirical output measures, particularly the mean absolute deviation of the relative errors of different forecasting methods. For each simulation run, we compared the forecasting of returns, net demands, their variance estimates, the reorder points, and the cost effects of inaccurate forecasts for each approximate forecast method (Methods 1, 2, and 4) relative to the exact Method 3, the latter using the maximal information. In each comparison we calculated the three aforementioned (mean positive and negative deviations, mean absolute deviation and histogram of deviations) evaluation measures but, for simplicity, here we shall only illustrate a single measure of relative error, namely the mean absolute deviation, defined for i = 1,2,4 by MAD(i) = mean of

Forecast(i)- Forecast (3)

(28)

Forecast (3) where the index i refers to the forecasting Methods 1, 2, and 4 described earlier, and the base case of comparison is Method 3.

Journal of Operations Management

27

An improvement (reduction) in the MAD is usually achieved by moving from a particular forecasting method to a more sophisticated one. The improvement can also be measured in many different ways. For our comparisons we have chosen to define the improvement factor for moving from method i to method k by IMP(i to k) = MAD(i) - MAD(k)

(29)

MAD(i) that characterizes the relative improvement achieved in MAD. Empirical Results of the Simulation Runs

In the following we describe the most important empirical results of the simulation experiments. 1. Forecast of Lead Time Returns from Previous Issues

The lead time returns from previous issues is the part of the overall forecast (lead time demand minus total returns in the lead time) that is most sensitive to the realizations of the recent issues and returns. Thus, for this quantity we see the most pronounced effects of the different forecast methods as functions of the system parameters such as the return distribution, the lead time and the demand parameters. The empirical tendencies will be described for MAD(2) defined by (28), which characterizes the forecast with the use of only information on recent issues. a) The return distribution pl ..... Pn has the following influence on the forecast error n • MAD(2) is monotonic increasing with n for fixed P = ~ Pi, as is illustrated by i=l example in Table 3 where MAD(2) is expressed in percentage form. • MAD(2) is monotonic increasing with _P for fixed n. • The shape of the distribution has no systematic effect on MAD(2). b) MAD(2) is monotonic increasing with the lead time L for L < n - 1, but independent of L for L > n - 1. (The latter result can also be proved analytically.) c) MAD(2) slightly increases with the coefficient of variation ~//x of demand. For the numerous simulation runs over a wide set of parameter values, the ranges of MAD(2) and MAD(4) were found to be Method 2 (just recent issues): 10-38% Method 4 (recent issues and aggregate returns): 4-20%. The one exception is the case where the return probability P = 0.4, but such a low value is unlikely to occur in practice. The average behavior Method 1 does not forecast the lead time return of previous issues separately, only the total return, thus MAD(1) could not be evaluated for lead time returns. For a given set of parameters there is a clear decreasing error tendency in moving from Method 2 to Method 4. This tendency is characterized by the improvement factor defined in (29), which took on values in the range: IMP(2 to 4): 19-71% which means a considerable relative improvement associated with moving from Method 2 to Method 4. The improvement factor slightly increases with n for given P and similarly slightly increases with P for given n. We have not recognized any systematic dependence of IMP(2 to 4) on other parameters.

28

TABLE 3 T H E E F F E C T S O F P A N D N O N MAD(2) (IN P E R C E N T ) O F L E A D T I M E R E T U R N S F R O M P R E V I O U S I S S U E S ( F O R L = n AND o'/tt=0.1) I'l

p

2

3

4

5

9

0.9 0.85 0.75 0.4

16.2 18.0 14.1 1.7

19.3 18. I 16.2 1.7

22 18.0 16.7 2.2

25 20 19 2.4

32 25 21 2.9

2. Forecast of Total Lead Time Returns For values o f L < n- 1 the same tendency to increase with L can be observed for the relative error of estimate of total lead time returns as was the case for the forecast of returns from previous issues. However, the increase is not as pronounced. Furthermore, for L > n-1 the relative error, MAD(2) defined by (28), is decreasing and tends to 0 as L goes to infinity. This latter property can be justified as follows. As we saw earlier, only the part of the total lead time returns that is from previous issues can be forecasted more accurately, using recent return data. This part remains unchanged for L > n - l , while the returns of lead time issues increase indefinitely with L, eventually completely dominating the result. The dependence on L is illustrated in Table 4. The tendencies for the other parameters of the simulations are the same as for the forecast of lead time returns from previous issues, detailed in the preceding subsection. The improvement factors were observed to lie in the following ranges IMP(I to 2): 25-76% IMP(2 to 4): 19-71% The latter range is seen to be the same as that found for the forecasts of returns from only the previous issues (in the preceding subsection). The equality of these two ranges can be proved analytically. TABLE 4 EFFECT OF LEAD TIME L ON THE RELATIVE ERROR (IN P E R C E N T ) O F F O R E C A S T O F T O T A L L E A D T I M E R E T U R N S n

I

2

3

4

2 3 4 5 9

4.8 5.2 7.2 6.6

5.6 6.1 9.1 9.5 18.2

5.5 10.4 14.8 15.3

4.9 5.7 12.7 15.5 18.7

L 5

8

I0

10.2

2.1 2.5 8.6 10.2 20.1

15.1 20.1

20

30

18.5

0.4 0.8 1.4 4.8 7.6

3. Reorder Point Estimation The reorder point is calculated as s = ED + k* X / V D

(5)

using the estimated values of lead time net demand ED and its variance VD, which are different for the four forecasting methods. The safety factor k* is an input parameter of the simulation experiments. The specific values indicated in Table 2 were used for k*.

Journal of Operations Management

29

For the relative error of reorder point estimation we observed the same tendencies in terms of dependence on the parameters of the return distribution and the demand as was the case for the forecast error of lead time returns from previous issues, detailed under subsection 1. We illustrate a part of these tendencies in Table 5 for MAD(2). In terms of the lead time, there is a tendency for an increase in relative error for L < n - l , then a decreasing tendency for L > n-1. However, these tendencies are not as pronounced as was the case for the relative error of forecasts of lead time returns, detailed in an earlier subsection. An illustration of the effect of L on the relative error MAD(2) is shown in Table 6. The observed ranges of the relative errors of different forecasting methods were as follows: Method 1 (average behavior): 10-250% Method 2 (just recent issues): 3-170% Method 4 (recent issues and aggregate returns): 0.5-30% For a given set of parameters errors are inclined to decrease as we move to the more advanced methods. This tendency is also reflected in the improvement factors which are in the ranges IMP(1 to 2): 15-98% IMP(2 to 4): 18-76% This means that a considerable improvement occurs in most cases in moving from Method 1 to Method 2 or from Method 2 to Method 4. The improvement factor slightly increases with P for given n and with n for given P No systematic dependence on the other system parameters was recognized.

TABLE 5 E F F E C T O F R E T U R N D I S T R I B U T I O N P A R A M E T E R S ON T H E R E L A T I V E E R R O R (IN P E R C E N T ) O F R E O R D E R P O I N T E S T I M A T I O N ( F O R L = 2) fl

p

2

3

4

5

9

0.9 0.85 0.75 0.4

26.3 22.1 20.0 2.0

48.3 36. I 42.1 1.2

52.1 40.2 36.4 1.7

40.6 41.4 38.2 3.5

112 71.2 66.4 4.8

TABLE 6 EFFECT OF LEAD TIME L ON THE RELATIVE ERROR (IN P E R C E N T ) O F R E O R D E R P O I N T E S T I M A T I O N ( F O R P = 0.9) L n

1

2

3

4

2 3 4 5 9

8.8 14.5

26.3 48 52 40.6 112

15 20 51.2

17.2

30

146

46 55

5

8

10

20

30

4.4 7.2 3.0

5. I

3.2 57 24 89

9.8 37

14

4. Cost Penalty of Using the Wrong Forecast The reorder point, as calculated in the form s = ED + k* V'-V-D ensures an optimal policy if the lead time net demand ED and its variance VD are estimated correctly. If these estimates are incorrect, the reorder point will be set too high or too low. The first alternative increases the inventory holding costs. However, the expected shortage costs decrease, but this decrease is less than the holding cost increase. Similarly, if the reorder point is set too low, the increase in the expected shortage cost is higher than the decrease in the expected holding cost. These properties are consequences of the choice of the appropriate safety factor k* which minimizes the expected total of the two costs. The increase in the expected total cost was calculated for each simulated realization of the past issues and returns for each forecast (mean and standard deviation) given by each of Methods 1, 2 and 4 relative to the base case estimation of Method 3. The means of the cost penalties serve as appropriate measures of the expected relative cost benefits of each of the more advanced forecasting methods. The mean cost penalty was examined as a function of the different system parameters. The parameters of the return distribution p~ ..... p, had the following effects: with an increase of n the mean cost penalty is increasing for fixed P = 2p,. However, with increasing P the cost penalty is decreasing for fixed n. The cause of this latter result is the low value of the lead time net demand when the return probability is very high (i.e., the loss probability is very low). In Table 7, we illustrate this property together with the decreasing tendency of the mean cost penalty with an increase in the lead time L. The value of the safety factor k* had no significant influence on the cost penalty for the relatively high service levels used in the experiments (see Table 2). The observed ranges of the mean cost penalties of the different forecasting methods were Method 1 (average behavior): 10-95% Method 2 (just recent issues): 2.5-52% Method 4 (recent issues and aggregate returns): 1-25% For any given set of system parameters a clear decreasing tendency was recognized in the cost penalty by moving to more advanced forecasting methods. This is demonstrated by the improvement factors being in the ranges IMP(I to 2): 30-96% IMP(2 to 4): 25-85% We did not observe any systematic dependence of the improvement factors on the system parameters. Finally Table 8 shows the average behavior of the different forecasting methods on a specific example using data obtained from industry. (These data on soft drink bottle issues and returns were published by Goh and Varaprasad (1986)). As an illustration, for L = 4, the MAD(l) for the reorder point estimation is 72.2 and IMP(I to 2) is 35.0. TABLE 7 E F F E C T S OF TOTAL R E T U R N P R O B A B I L I T Y AND LEAD T I M E ON T H E MEAN COST PENALTY (IN PERCENT) OF USING M E T H O D 2 INSTEAD OF M E T H O D 3 (FOR n = 5) L P 0.95 0,9 0.6

2 22.1 51.3 65.9

Journal of Operations Management

5 20.4 44.7 48.9

I0 15.6 34.6 45.4

30 10.1 23.9 29.1

31

RELATIVE ERROR

TABLE 8 AND IMPROVEMENT FACTOR (BOTH IN PERCENTS) OF THE DIFFERENT FORECASTS

Return distribution: Pl = 0.639, P2 = 0.084, P3 0.238 Coefficient of variation of demand: ~r/~ = 0.22 =

FORECAST

L =

2

4

8

MAD(I) MAD(2) MAD(4)

IMP(2 to 4)

no forecast 17.2 55.5 7.8

no forecast 17.2 55.5 7.8

no forecast 17.2 55.7 7.8

2. Total lead time returns

MAD(I) MAD(2) MAD(4)

IMP(I to 2) IMP(2 to 4) IMP( 1 to 4)

9.9 4.9 2.2

50.3 55.5 77.9

5.6 2.5 I. 1

54.6 55.6 79.8

2.8 1.3 0.6

54.6 55.7 79.9

3. Reorder point

MAD(I) MAD(2) MAD(4)

IMP(I to 2) IMP(2 to 4) IMP(I to4)

74.9 63. I 21.4

15.8 66.0 71.4

72.2 46.9 15.0

35.0 68.0 79.2

79.9 37.2 13.4

53.5 64.0 83.2

4. Cost penalty

MAD( I ) MAD(2) MAD(4)

IMP( 1 to 2) IMP(2 to 4) IMP( 1 to 4)

33.6 11.3 2.7

88.3 76.2 97.2

97.6 10.5 2.5

89.3 76.2 97.4

95.9 9.8 2.3

89.7 76.3 97.6

1. Lead time returns from previous ISSUES

CONCLUSIONS T h e results on a w i d e range of simulated data, including s o m e cases based on empirical data obtained from industry, s h o w e d that the relative errors in the estimated values are considerable w h e n one uses the simplest forecasting m e t h o d (Method 1). The cost penalties of using these incorrect forecasts are also quite high. The i m p r o v e d forecasting methods d e v e l o p e d in this paper, despite the mathematical c o m p l e x i t i e s in their d e v e l o p m e n t , are in fact quite easy to use and do not involve any c o m p l e x c o m p u t e r routines. The i m p r o v e m e n t s a c h i e v e d by applying each m o r e advanced m e t h o d are substantial. The data used by two o f the m o r e a d v a n c e d m e t h o d s ( M e t h o d 2 based on recent issues and M e t h o d 4 based on recent issues and aggregate returns) should be obtainable without substantial cost in most industrial settings. In contrast, the b e n c h m a r k best case m e t h o d (Method 3 based on recent issues and individually identified returns) involves detailed information on individual returns that is likely not obtainable without substantial cost, only justified for very e x p e n s i v e containers. ACKNOWLEDGEMENT The research leading to this paper was supported by the Natural Sciences and Engineering Research Council of Canada under Grant A1485. The authors would like to express their thanks to the editor and an anonymous referee for providing a number of constructive suggestions.

REFERENCES I. Flies, R.H., Jr., and Sr. Problems in Computing Trippage. Beverage WorM, vol. 12, 1977, 146-149. 2. Fujita, S. Application of Trippage Distribution to Recycling Problems. Proceedings of 5th International Symposium on Management Science, Osaka, 1977, 1-12. 3. Goh, T.N., and N. Varaprasad. A Statistical Methodology for the Analysis of the Life-Cycle of Reusable Containers. liE Transactions, vol. 18, 1986, 42-47.

32

4. Hayya, J., D. Armstrong, and N. Gressis. A Note on the Ratio of Two Normally Distributed Variables. Management Science, vol. 21, 1975, 1338-1341. 5. Hmkley, D.V. On the Ratio of Two Correlated Normal Random Variables. Biometrica, vol. 56, 1969, 635-639. 6. Marsaglia, G. Ratios of Normal Variables and Ratios of Sums of Uniform Variables. Journal of the American Statistical Association, vol. 60, 1965, 193-204. 7. Silver, E.A., and R. Peterson. Decision Systems for Inventor), Management and Production Planning, 2nd ed. NY: John Wiley & Sons, 1985. 8. Weymes, E. A Report on Trippage Rates of Carbonated Soft Drinks. Cranfield, U.K.: The Retail Development Centre, Cranfield School of Management, 1978.

APPENDIX A THE MIXED BINOMIAL M O D E L P r o p o s i t i o n 1: For a mixed b i n o m i a l r a n d o m variable b with a random n u m b e r n of trials and k n o w n success probability p E(b) = p E(n)

(A.l)

and Var(b) = p2 Mar(n) + p ( 1 - p ) E ( n ) .

(A.2)

Proof: Each trial Xm has a probability of success p, independently of one another for m = 1 ..... n. Thus P(Xm = 1) = p, P ( X m = 0 ) = q = 1-p and the generating function of Xm is given by f(s) = q + ps. co

Let P(n = 1) = gi and g(s) =

X

g, s' be the generating function of the r a n d o m variable n.

i=l The mixed b i n o m i a l b has the distribution: hk = P ( b = k )

=

P ( n = i ) P(X~ + X2 + ... + Xi = k), i=O

its generating function is: OO

X g, fi(s) = g(f(s)) h(s) = i=O and the expected value is: E(b) = h'(1) = p g ' ( q + p )

= p g'(1) = p E(n).

Further, the variance Var(b) = h"(l) + h ' ( 1 ) - [h'(1)] 2 = h"(l) + p E ( n ) - p2 E2(n) where h"(1) = d

[p g ' ( q + p s ) ] s = ~

= p2 g , , ( q + p ) = p2 g"(1).

ds

Since Var(n) = g" (1) + g ' ( l ) - [q' (1)] 2 = g"(1) + E(n) - E2(n),

Journal of Operations Management

33

it follows that h"(1) = p2 Var(n) - p2 E(n) + p2 E2(n) and thus Var(b) = p2 Var(n) + (p_p2)E(n) = p2 Var(n) + p(1-p)E(n).. Proposition 2: For a mixed binomial random variable b with random n and known p Var(n-b) = (l-p) 2 Var(n) + p(1-p)E(n)

(A.3)

Proof: For gi = P(n = i), i

E(nb) =

=

oo E i=l

i g~

00

.~. k=l

i E k (~) k=l

Y-. i gi k (ik) pk(1-p) '-k i=l pk(l_p),-k =

oo E i 2 gi P = P E(n 2) i=l

= p[Var(n) + E2(n)] thus Cov(n,b) = E(nb) - E(n)E(b) = p Var(n) + p E2(n) - p E2(n) = p Var(n) fu~her Var(n-b) = Var(n) + Var(b)- 2 Cov(n,b) = (l-p) 2 Var(n) + p ( l - p ) E ( n ) . .

APPENDIX B THE MULTINOMIAL MODEL

Consider a knownn return distribution p. ..... p, of each particular unit issued, with a probability p~ = 1 - j ~ l p, of never returning and with the assumption that the probability of the return of a unit is independent of the return of other units. Then the number of containers, v.j, issued in period i and returned in period i + j (j = 1..... n) is a sequence of multinomial trials with the multinomial probability distribution P(vi., = k, ..... v , . . =

n k..v,~ = ui-Z j=l

kj)

n u , - ~ kj k, k. j= 1 p, ...p. p~ ,

u,! = n

k, ! . . . k . ! ( u,- Y~ kj)! j=l where vi~ denotes the number of containers that never return from the issues ui of period i. The total return W, in the lead time periods t + 1..... t + L from issue u, is the sum of random returns in each of these periods: Wl

34

=

Vt,t-i + 1 +

Vt.l-I + 2 "[- ' ' '

"1- VI,t-I + L

for i < t

and W i = Vi,l + vi,2 +

for t < i < t + L.

... + viA-*+ L

Since the vi.j (j = 1..... n) have a multinomial distribution for each i > 1, thus Wi, being the sum of multinomial components, has a binomial distribution where u~ is the number of trials and the success probability is the sum of the appropriate component probabilities (pj's). We denote each success probability by Ri(pj's). Then, because the component probabilities, the pj's, are 0 for j < 1 and for j > n, R~(pj's) = 0 jm ~

Ri(pj's) =

for i < t-n + 1, Pt-i +j

for i m < i < t,

(B.l)

j=l with i m = andR~(pj's) =

max{1,t-n+ 1} and jm = min{L,n+i-t}, jn Z pj j=l

fort
with jn = min{n,t-i+ L}. The total number of containers returned in period m m-I Ym =

]~

Vi,m-i

i=m-n is the sum of returns from independent issues, however, the returns in the different periods are correlated. We next develop the components of the covariance matrix (T) of the y's. For the multinomial distribution COV(Vi.j,Vi,k) = - Ui pj Pk, t-k thus Tj.k = C o v ( y t _ j + l , Y t _ k + l ) = -

~ Ui Pt-l+l-i Pt-k+l-,

i = ij

for j = 1..... n-2, j < k < n-2 and ij = max{l,t-j + l-n},

(B.2)

also Tk.j = Tj.k for k < j

t-j further T~ = Var(yt_j + l) =

i=ij

ui Pc-j+ J-i (1 - Pt-j + ~-,)

for j = 1 ..... n-2. The total lead time return from previous issues: t-1 E W, im

with im = m a x { l , t - n + 1}

i=

is correlated with Yt-j+] for j = 1 ..... n-2 and according to the multinomial model (as above) we find

(t~l Cj ~ C o y

) Wi, Yt-j+l

n-j = -

".i =im

E

i=im

jr, I.lt_n+i Pn-j+l-i

]~

Pn-i+m

(B.3)

m= 1

for j = I ..... n-2 and jr, = min{i,L}.

Journal of Operations Management

35