Robust compound Poisson parameter estimation for inventory control (cid:2)

Most companies store demand data periodically and make periodic demand forecasts, whereas many de- mand processes in inventory control need parameter estimates at the individual customer level. Guidance on estimating the parameters of a continuous-time demand process from period demand data is lacking, in particular for the popular and well-studied compound Poisson class of demand. Whereas the statistics literature typically focuses on asymptotic properties, parameters for inventory control have to be estimated based on a limited number of periodic historical demand observations. We show that the standard Method-of-Moments (MM) estimator – the default choice in applied inventory control research – is severely biased for ﬁnite samples. The Maximum Likelihood (ML) estimator – which needs to be obtained by a numerical search – performs better, but both estimators lack robustness to misspeciﬁcation of the demand size distribution. We propose an intuitive, consistent, closed-form MM alternative that dom-inates standard MM and ML in terms of estimation accuracy and on-target inventory performance. Its closed form does not depend on the speciﬁc demand size distribution, making it robust and easily appli-cable in large-scale applications with many items. In a case study, we ﬁnd that the accuracy loss due to storing demand periodically is four times as high under standard MM as under the proposed estimator.


Introduction
The compound Poisson demand process is one of the most general and widely studied demand processes for inventory control. It has a natural business interpretation with customers arriving in continuous time, and their demand sizes following any general compounding distribution, and is therefore a standard choice in textbooks (e.g. [1][2][3] ) and commercial software (e.g. Forecast Pro, Oracle NetSuite, SAP APO, SAS, and Slimstock), especially for intermittent demand patterns. After Galliher et al. [4] studied the special case of Poisson-geometric demand, Feeney and Sherbrooke [5] and Chen et al. [6] derived order-up-to policies under general compound Poisson demand. Archibald and Silver [7] extended this work to general (s, S) policies. Order-up-to policies can be deter-mined efficiently under compound Poisson demand, as shown by Sherbrooke [8] and Graves [9] , who provided approximate procedures, and Axsäter [10] , who presented an exact method. Several empirical studies concluded that, especially if products are not fast-moving, compound Poisson distributions provide the best fit (e.g. [11][12][13] ). Therefore, they remain a standard choice in inventory control research (e.g. [14][15][16] ).
In order to exactly characterize inventory policies, a complete specification of the compound Poisson demand process is required. This is most obvious for continuous-review inventory systems. Here, the arrival of a customer immediately triggers an order, and the probability distribution of individual customer demand sizes determines the expected reorder level undershoot. However, also under periodic review all parameters of the compound Poisson demand process are needed for exact calculations [17] . As compound Poisson demand processes operate in continuous time, it would be natural to estimate their parameters from individual transaction data. The average number of customer arrivals in a period is an unbiased estimator for the arrival rate and the parameter(s) of the compounding distribution can be estimated directly from the demand sizes. However, it is common in practice (and a typical assumption in applied inventory control literature [12,13,18]  only total demand per period is stored and used as a basis for forecasting. In fact, temporal data aggregation to improve accuracy is frequently suggested (see e.g. [19][20][21] ). As a result, commercial forecasting software is mainly based on period demand. Storing demand per period rather than per individual transaction does lead to an accuracy loss. As we show in this paper, both the choice of the parameter estimator and the period length to which demands are aggregated have a large influence on the magnitude of the error.
Separate forecasting of demand sizes and demand inter-arrival times may be achieved by so-called Size-Interval methods. The most well-known one -twice applying exponential smoothingwas first proposed by Croston [22] . A survey of software packages' forecasting capabilities [23] confirms the popularity of Croston's method. Continuations of this research stream are reviewed by Syntetos et al. [24] . However, period forecasting methods estimate the mean demand size in a period and the mean interval between periods with demand occurrence, rather than the mean demand size of an individual customer and mean interval between individual customer arrivals. Therefore, they cannot be directly used in inventory control, especially not if demand is intermittent [25,26] . Ignoring this mismatch leads even asymptotically to significantly biased estimates and dramatically underachieved service levels, as we will show in Section 2 .
Very few authors in the inventory control literature discussed how the unknown parameters of compound Poisson demand processes (or of demand processes in general) can be obtained from period demand data. The authors who did, suggested the use of standard Method-of-Moments (MM) estimators [2,27] , which are therefore at present the standard choice in applied work [13,18] . In the wider (statistics) literature, the discussion of parameter estimation centers around (variants of) MM and Maximum Likelihood (ML). Examples are Anscombe [28] and Sichel [29] (logarithmic compounding distribution), Öztürk [30] (exponential compounding distribution), Withers and Nadarajah [31] (gamma compounding distribution), and Balakrishnan et al. [32] (geometric compounding distribution). MM and ML variants are also given for certain compound Poisson distributions in textbooks such as Johnson et al. [33] , without a performance comparison. If authors compare properties of different estimators, then the discussion is almost always limited to large-sample properties. In this paper we focus on finitesample performance instead, which is important in practice. Product life-cycles are shortening, and even if historic demand data are available for multiple years, then not all of them may be indicative of future demand.
Our main contribution is the presentation of an alternative MM estimator that retains Croston's core idea of separately estimating customer arrival rates and demand size means. It bases the customer arrival rate on the fraction of periods without demand, and then uses central moments to estimate the parameters of the (individual) demand size distribution. This estimator has an intuitive, closed-form solution, independent of the specific compounding distribution. It is furthermore robust to misspecification of the compounding distribution, which is an important practical advantage. The idea of estimating distribution parameters based on the 'zero frequency' has sporadically been applied [12,28,[34][35][36] . However, i) a general method to estimate any compound Poisson demand parameters based on this logic has never been given, and ii) a study of such a method's performance is lacking.
Focusing on the exponential compounding distribution in the main body of the paper and discussing the geometric compounding distribution in Appendix B , we compare the finite-sample bias and variance of standard MM, ML, and the proposed alternative MM estimator. Thereafter, we compare the achieved fill rates in a continuous review order-up-to system. We furthermore study all estimators under misspecification of the demand size distribution and find that both standard MM and ML lack robustness, whereas the proposed method is unaffected by the specific demand size distribution. A case study validates the results and quantifies the accuracy loss due to storing demand per period rather than per transaction.
The remainder of this paper is organized as follows. Section 2 introduces the inventory model, demand model, fill rate service measure, and discusses the misapplication of period Size-Interval methods. Section 3 introduces the consistent standard MM and ML estimator and the proposed estimator. Section 4 provides a comparison of the consistent estimators' biases, whereas Section 5 discusses their variances. Section 6 presents a comparison of achieved fill rates in an order-up-to inventory model, and Section 7 gives robustness results. A comparison of the estimators in a real case study and a quantification of the accuracy loss due to storing aggregated demands per period rather than individual transactions is given in Section 8 . Section 9 concludes. Appendix A contains derivations related to Section 4 , and Appendix B presents results for the geometric compounding distribution.

Demand model, inventory model, and size-interval methods
In this section we specify the general compound Poisson demand model and the compounding distributions on which we base our results. Thereafter we discuss the order-up-to inventory model and the fill rate calculation. Finally, we discuss the consequences of misusing period Size-Interval estimates as compound Poisson demand parameters.

Demand model
We consider a model where demand follows a compound Poisson process with some general (discrete or continuous) compounding distribution. We have a sample of n period demand realizations x 1 , ..., x n . That is, for t = 1 , 2 , ..., n , we observe realizations of the random variable X t = K i =1 D i , where K follows a Poisson distribution with mean λ and, independently, the individual customer demands D i are i.i.d. with mean μ c . The individual customer demands are not observed. We require that P (D i ≤ 0) = 0 as a customer can naturally only have positive demand. The set x 1 , ..., x n of period demand realizations is used to estimate the parameters of the compound Poisson process.
In the remainder of the paper we focus on the exponential compounding distribution which is continuous and therefore more suitable for demonstration, since any fill rate level can be achieved exactly. We refer to Appendix B for additional results based on the geometric compounding distribution. For the exponential distribution, the individual demand sizes D i have the probability density function for x ≥ 0 and θ > 0 , with demand size mean μ c = 1 /θ . For the geometric distribution, the individual demand sizes D i have the (discrete) probability mass function for x = 1 , 2 , ... and 0 < β < 1 , with demand size mean μ c = 1 / (1 − β ) .

Inventory model
We study a continuous review order-up-to inventory control model with a positive fixed lead time of L periods and full backordering, under a fill rate constraint. Order-up-to inventory models are relatively easy to optimize, popular in practice, and theoretically optimal for various cost and service models when order costs are negligible [1,37] . As a service target we use the fill rate, which is the expected fraction of a demand that can be serviced immediately from on-hand inventory.
We follow the computation of the fill rate for compound Poisson demand models by Axsäter [2] , and generalize it to handle also continuous compounding distributions. Specifically, denoting an individual demand by D i and the inventory level by IL , for a given order-up-to level S, the fill rate is the expected fraction of a single order that can be fulfilled from on-hand stock, The distribution of D i is the selected compounding distribution with mean μ c , whereas the distribution of IL follows from the observation that the inventory level at time t + L is equal to S minus the lead time demand. Since demand per period is compound Poisson with rate λ and demand size mean μ c , the demand during the lead time follows a compound Poisson distribution with the same compounding distribution and rate λL . Denoting lead time demand by D L , we find the relationship P ( where γ is the target fill rate. In the remainder of this paper, we use γ = 0 . 95 (i.e. 95%).

Misapplication of period size-interval methods
Period Size-Interval methods yield estimates of the time between two periods with positive demand and the size of a nonzero period demand. We demonstrate here the consequences of misapplying these as compound Poisson parameter estimators. We first discuss Croston's method and a 'bias-corrected' version of this method. Thereafter we discuss an intuitive method that uses the average of all non-zero period demands and the average of all times between periods with positive demand, which we refer to as 'Unweighted Averaging'. All discussed methods only update after a period with positive demand.
Croston's method estimates in period n the time between two periods with positive demand (denoted by ˆ p n ) and the size of a non-zero period demand (denoted by ˆ s n ) both by exponential smoothing. Formally, given that we have observed n periods, the n th period has a positive demand x n , and the previous positive demand was q n periods ago, where 0 < α 1 < 1 and 0 < α 2 < 1 are constants.
If Croston's period estimates are wrongly taken as the estimates of the individual demand parameters, then ˆ s n is the estimate after n period observations for the mean μ c of the compounding distribution, whereas 1 / ˆ p n is the estimate for λ. As shown by Syntetos and Boylan [38] , using 1 / ˆ p n as an estimate for the arrival rate of a compound Poisson process leads to an inversion bias, as E [1 / ˆ . This inversion bias does not diminish as the sample size goes to infinity, as Croston's method uses exponential smoothing and thus does not use all data points, so that the Central Limit Theorem does not hold. Syntetos and Boylan [39] approximately corrected for this bias using a second-order Taylor expansion, and derived the estimator (1 − α 2 / 2) / ˆ p n for λ. The esti-mator for μ c does not suffer from the inversion bias and remains unaltered.
It is straightforward to show that, with or without the inversion bias correction, Croston's method yields inconsistent customer demand parameter estimates when applied to period demand. Consider ˆ s n , which is unaffected by the inversion bias. Instead of estimating the size of an individual demand, it estimates the size of a positive period demand. Therefore, , which is strictly larger than μ c and smaller than μ c λ for all λ > 0 and μ c > 0 . This limit approaches μ c as λ → 0 , and μ c λ as λ → ∞ . The estimator's limit (and thus its bias) is proportional to μ c .
The asymptotic bias of the arrival rate estimate is more cumbersome to quantify. Instead of estimating the time between two customer arrivals, ˆ p n estimates the number of periods between two periods with positive demand. This number of periods (denoted by q ) follows a geometric distribution with success probability ) as n → ∞ , we can only approximate the limit of E [1 / ˆ p n ] . We do so by a second-order Taylor expansion and find, using the mean and variance of the geometric distribution and the fact that ˆ p n is an exponential smoothing estimator with asymptotic variance Var As α 2 → 0 , this limit approaches 1 − exp (−λ) , which is the result without inversion bias. In that scenario, λ is strictly underestimated for all λ > 0 . The approximate limit of the estimator converges to 0 as λ → 0 , and to 1 as λ → ∞ . Obviously, the bias correction proposed by Syntetos and Boylan [39] scales the estimated arrival rate, and thus its expectation, by (1 − α 2 / 2) . In conclusion, Croston's method generally underestimates the arrival rate, but this is somewhat counteracted by its inversion bias, leading to an overestimation especially for low values of λ and/or high values of α 2 .
However, (approximately) correcting for that inversion bias intensifies the underestimation. While Croston's method and its variants give greater weight to the most recent demand sizes and intervals, this is not a requirement for Size-Interval methods. An alternative approach is to use equal weights, but retain the separate estimation of demand size means and intervals. The Unweighted Averaging method does this by taking as the demand size mean the average of all positive period demand sizes, and as arrival rate the inverse of the average number of periods between two periods with positive demands. Given that we have observed n periods, of which n p had a positive demand, the n −th period has a positive demand x n , and we have observed the previous positive demand q n periods ago, Unweighted Averaging sets It is straightforward to observe that, just as with both Croston variants, . Unlike Croston's method, the Unweighted Averaging method is asymptotically unbiased, since it uses all positive demands and inter-arrival times in the data set rather than applying diminishing weights via exponential smoothing. Therefore, the Central Limit Theorem holds and ˆ p n converges in probability to E [ ˆ p n ] , so that The Unweighted Averaging estimator for λ converges to 0 as λ → 0 , and to 1 as λ → ∞ . This reveals the flaw in using period demand estimators as individual demand parameters: no distinction can be made between 1 or several demands in a period, leading to an underestimated arrival rate. All individual demands in a period are added and considered as a single customer demand, leading to an overestimated demand size mean. Table 1 quantifies the (approximate) asymptotic bias resulting from using standard Croston's method ('CROST'), the corrected version by Syntetos and Boylan [39] (Syntetos-Boylan Approximation, 'SBA'), or the Unweighted Averaging method ('UA'), for several parameter choices. Furthermore, we report the fill rate that would be achieved for both the exponential and geometric compounding distribution, under the true demand parameters, when a target fill rate of 95% is chosen. We use a lead time of 2 periods.
None of the results depend on μ c , except for the achieved fill rates under the geometric compounding distribution. The reason is that for discrete demand distributions, fill rates cannot be achieved exactly. For low arrival rates and/or high smoothing parameter values, Croston's inversion bias results in an overall overestimation of λ. The Syntetos-Boylan Approximation achieves an asymptotic underestimation for all scenarios. The Unweighted Averaging method uniformly underestimates λ. The demand size mean μ c is always overestimated for all methods, and all (absolute) percentage asymptotic biases are increasing in λ.
Interestingly, all achieved fill rates are asymptotically too high, irrespective of whether λ is over-or underestimated. This is because the achieved fill rate depends more strongly on μ c , as we will further discuss in Section 6 . The achieved fill rates are increasing in λ, and for λ = 1 they are over 97.8% for all compounding distributions. We conclude that misapplying period Size-Interval methods to estimate compound Poisson demand parameters leads to severely overestimated demand size means, significantly overshot fill rates and correspondingly excessive inventory costs.

Consistent estimators
Having discussed the misapplication of period forecasting methods for obtaining compound Poisson demand parameters, we move to consistent estimation methods. We revisit the compound Poisson demand model introduced in the previous section and discuss the standard MM estimator (applied in the inventory control literature), the ML estimator, and the proposed alternative MM estimator -first for general compound Poisson demand processes and then for the specific case of an exponential compounding distribution.
We first give some properties of the Poisson-exponential distribution that are helpful for our analysis [30] . Its probability density function is given by for λ > 0 , μ c > 0 , and x > 0 . This distribution has a positive probability mass of exp (−λ) at x = 0 . The first four central moments are μ 1 = λμ c , μ 2 = 2 λμ 2 c , μ 3 = 6 λμ 3 c , and μ 4 = (24 λ + 12 λ 2 ) μ 4 c . We assume that a sample x 1 , ..., x n of period demands is available and denote the i th central sample moment of the period demand observations by ˆ μ i . In particular, ˆ μ 1 denotes the sample mean and ˆ μ 2 the sample variance.

The standard method-of-moments estimator
Standard MM estimators equate central moments of the compound Poisson demand distribution to their corresponding sample equivalents. That is, if the compounding distribution has P parameters, we solve for the unknown arrival rate and the P parameters of the compounding distribution. An exponential compounding distribution has one parameter, which we take to be its mean. So, there are in total two parameters to be estimated. The standard MM estimators for the arrival rate and the demand size mean are then MM estimators are guaranteed to be consistent, but may still have severe finite-sample biases and variances, as our numerical results will reveal. Furthermore, the second order sample moment (the variance of the compound Poisson distribution) depends on the exact shape of the compounding distribution and not only on its mean. This implies that standard MM is not robust to misspecified demand size distributions.

The maximum likelihood estimator
The ML estimator of a compound Poisson distribution is found by maximizing the log-likelihood function with respect to the unknown parameters. In most cases, the maximizing arguments cannot be found in closed form, and a numerical search procedure has to be invoked. The log-likelihood function is defined as where f CP denotes the probability density function of the corresponding (full) compound Poisson distribution and θ denotes its vector of parameters, including the arrival rate of the Poisson process and the parameters of the compounding distribution.
For the exponential distribution, we have f CP = f PE and θ = (λ, μ c ) . One can maximize L (λ, μ c | x 1 , ..., x n ) to find estimates ˆ λ and ˆ μ c , using a non-linear maximization algorithm. Although ML is computationally more demanding than the other methods, it is not prohibitively expensive, as it involves only a two-dimensional search. Moreover, it is guaranteed to provide a consistent estimator that asymptotically has the lowest variance. However, ML offers no guaranteed optimal performance for finite samples.

The proposed method
We propose an alternative MM estimator and choose its moments such that the intuition of the Size-Interval forecasting literature -separate estimation of the customer arrival rate and the demand size mean -is employed. Also this estimator is guaranteed to be consistent, as it is an MM estimator. We first discuss estimation of the arrival rate λ. Denote by N 0 the (stochastic) number of periods in the sample of size n in which no demand occurred. The corresponding realization of that number of periods is denoted by n 0 . The occurrence of no demand in a period is Bernoulli distributed with probability exp (−λ) , and therefore the number of periods without demand out of n total periods follows a binomial distribution with parameters n and exp (−λ) . Therefore We equate this sample moment to its expectation to find the estimator for the arrival rate: For each of the parameters of the compounding distribution, we equate (as with standard MM) central moments of the compound Poisson distribution to their sample equivalents. However, the explicit use of the fraction of periods without demand to estimate the arrival rate implies that the highest order of moments used by this method is always one lower than for the standard MM approach. That is, if the compounding distribution has P unknown parameters, we solve ˆ λ = − ln (n 0 /n ) , In the case of an exponential compounding distribution (more generally, in the case of any compounding distribution with a single parameter), the set of equations reduces to ˆ λ = − ln (n 0 /n ) , μ 1 = ˆ μ 1 .
Using the fact that μ 1 = λμ c , we solve for the demand size To calculate these estimates, one only needs to obtain the fraction of periods without demand and the overall mean demand during all periods. Furthermore, whereas standard MM and ML have to be derived separately for every compounding distribution, this method gives estimates that do not depend on the specific compounding distribution. Irrespective of the compounding distribution, ˆ λ and ˆ μ c have exactly the same closed-form solution. In Section 7 we elaborate on this property, which extends to compounding distributions with an arbitrary number of parameters. Furthermore, we show that this property implies robustness: irrespective of the compounding distribution, the proposed estimators converge to the corresponding true parameters.

Bias comparison
In this section we compare the biases of the three estimators discussed in the previous section. Asymptotically, these estimators are all unbiased, but for finite sample sizes this is not the case. As standard MM and the proposed estimator both exist in closed form, we can analytically study their biases, as we will do via a Taylor approximation, and find the conditions under which the proposed estimator is more accurate. Thereafter we perform a numerical study to also compare both estimators with ML.

Analytical bias comparison of standard MM and the proposed estimator
To analytically compare standard MM and the proposed estimator, we approximate the estimators by a second-order Taylor series, as a function of the sample mean ( ˆ μ 1 ) and the sample variance ( ˆ μ 2 ), and then calculate the expected value. For the derivation, we refer to Appendix A . Here we study the results, starting with the arrival rate estimator. For standard MM, we find: For the proposed estimator, we find: For the proposed estimator, we find: These estimators both have a negative bias for all values of λ, μ c , and n . The demand size mean μ c acts solely as a scaling parameter. It follows that for λ < 2 the new estimator is more accurate than standard MM. The scenario with λ = 2 corresponds to a 14% probability of zero demand in a period. We see similar dynamics as for the arrival rate estimator. For λ = 1 the proposed estimator is approximately 7 times as accurate as standard MM, whereas for λ = 1 / 4 this factor has increased to 170.
There are many real life settings with moderate to high intermittency. Lengu et al. [12] report that more than 80% of their studied items have an estimated arrival rate below 1, the majority of which are below 1/4. Turrini and Meissner [13] report median demand intervals of 10 and 17 periods, and a maximum of 80 periods. In Section 4.2 we present more insights into the biases of standard MM, ML, and the proposed method by means of a numerical analysis.

Numerical study of all estimators' biases
In this section we perform a numerical analysis of the biases of standard MM, ML, and the proposed method, under an exponential compounding distribution. We select λ = 1 / 16 , 1 / 4 , and 1, representing high, moderate, and light intermittency, respectively, and present the results for μ c = 2 only. Experiments with other values and the analysis in Section 4.1 show that μ c acts only as a scaling parameter. We let the sample size range from n = 4 to n = 200 . Each time we draw 1,0 0 0,0 0 0 period demand histories from the corresponding compound Poisson distribution. For each draw, we calculate the parameter estimates for λ and μ c using (i) the standard MM estimator, (ii) the ML estimator, and (iii) the proposed method. We present per parameter combination, per sample size, and per estimator, the average of all obtained estimates.
Two scenarios deserve special attention. The first is the case where only periods without demand are observed. In practice one will not fit any demand distribution then, but in numerical experiments it does occasionally occur. All methods are then undefined and to enable a fair comparison, we set in this case all estimated arrival rate to 0 and do not define demand size mean estimates.
The other extreme scenario is that with no periods without demand. In this unlikely event, we set the proposed method equal to standard MM. Fig. 1 shows the results, with on the left hand side the estimated arrival rates λ and on the right hand side the estimated demand size means μ c . Note that in this and all following figures the marker locations have been shifted between estimators to enhance readability. The scenario with λ = 1 / 16 ( Fig. 1 a and b) corresponds to a 94% probability of having zero demand in a period. The standard MM estimator converges the most slowly by far. For n = 200 , its arrival rate estimate is still 13% too high, whereas its demand size mean estimate is 8% too low. The proposed method and the ML estimator perform almost identically, and are considerably more accurate than standard MM. For n ≥ 50 the difference with the true parameters is below 1%. As the sample size decreases, the probability of observing no demand at all increases.
Especially for low arrival rates, this affects the results.
The scenario with λ = 1 / 4 ( Fig. 1 c and d) corresponds to a 78% probability of having zero demand in a period. Standard MM still shows a significant deviation from the true value for n = 200 (5% for λ and 2% for μ c ). The arrival rates converge more slowly than the demand size means, and the proposed estimator actually (slightly) outperforms ML. In Fig. 1 e and f, the arrival rate is increased to 1, corresponding to a 37% probability that a period has zero demand. The standard MM estimator still performs much worse than the proposed method and ML, but the gap is smaller than for lower arrival rates. The proposed method now outperforms ML by a larger margin.
Comparing relative accuracies, we find that the ML and proposed arrival rate estimator perform best for low λ, whereas the standard MM arrival rate estimator becomes better when λ increases, although it is still severely outperformed even for λ = 1 . All estimators for the demand size mean are more accurate as λ increases, as fewer periods without demand occur. Directly related to this observation is the question of the optimal period length (e.g. hour, day, week, or month). Different period lengths lead to different levels of intermittency, corresponding to different true values of λ. In Section 6 we study the resulting effect on inventory per- formance. An explanation for the very similar performance of ML and the proposed estimator for low arrival rates (here and for the Poissongeometric distribution in Appendix B ) can be found in their definitions. The probability mass of a compound Poisson distribution for periods with no demand reduces to exp (−λ) , which is minimized at λ = 0 . For low arrival rates many periods have no demand. The sample log-likelihood function is the sum of the log-densities of all periods and therefore relatively insensitive to (few) nonzero observations. This insensitivity also holds for the proposed arrival rate estimator, as the natural logarithm is relatively flat around 1 (the point where n 0 = n ). The standard MM arrival rate estimator, however, is more sensitive to (even few) nonzero observations, as it uses the squared mean and the variance of the sample. For higher arrival rates, the demand pattern becomes less intermittent, benefiting standard MM. That is also where the bias of ML differs more from that of the proposed method (e.g. Fig. 1 e and f).
The results of this section can be used to validate the accuracy of the analytical bias approximation of Section 4.1 . Table 2 presents the comparison for the case λ = 1 / 4 , μ c = 2 . The proposed estimator's expectation is better approximated by the second-order Taylor series than the standard MM estimator, as would be expected based on their different functional forms. For n = 10 , the analytical approximation for standard MM differs substantially from the numerical evaluation. There are two reasons for this: 1) the Taylor series approximation is generally more accurate for larger samples, and 2) for small sample sizes the special cases mentioned above occur more frequently. For n ≥ 50 the difference between the an- alytical and numerical results is small. The proposed estimator is well approximated by the Taylor series for all sample sizes. Furthermore, under both the analytical and numerical evaluation, the methods retain the same ranking.
The proposed method has another important advantage over ML: it has a simple, closed-form solution, whereas optimization of the log-likelihood function for ML involves a two-dimensional numerical search, each time evaluating a function with infinite sums. Although this procedure is not prohibitively expensive, it does require much more computation time than standard MM or the proposed estimator. Table 3   takes 61 minutes by ML, and less than 3 seconds by the proposed or standard MM estimator.

Variance comparison
In this section we use the same numerical set-up as in the previous section, but report the variances over the draws. Fig. 2 shows the results, with the arrival rates on the left hand side and the demand size means on the right hand side. We firstly observe that the proposed method and ML perform very similarly. This can be explained analogously to their similar biases (see Section 4.2 ): they differ only in their treatment of non-zero observations, to which both are relatively insensitive compared to the standard MM estimator. Since the variance of ML acts as an asymptotic lower bound and it is practically indistinguishable from the proposed estimator, we conclude that the new estimator shows near-optimal performance.
For the arrival rate estimators, standard MM has a substantially larger variance than the other two methods. As the fraction of nonzero observations increases with λ, so does the difference between standard MM and the other two methods. However, for the demand size mean it can be observed that only from a certain minimum sample size onward, standard MM has a larger variance than the other two methods. This reflects a drawback of our numerical set-up. If only periods without customer arrivals are drawn, the demand size mean is not estimated. The reported averages are therefore biased towards 'less intermittent' draws, favoring standard MM. For larger sample sizes (where standard MM does have a significantly larger variance than the other methods) this event virtually never occurs.

Inventory control implications
For the order-up-to inventory model described in Section 2.2 , we study the service effects when demand parameters are obtained according to the three consistent estimators, again for an exponential compounding distribution. We use a target fill rate of 95% and report the achieved fill rates. Fig. 3 a shows the results for λ = 1 / 16 , μ c = 2 , and L = 2 . Consistently with the results of Section 4 , the ML and the proposed method converge much faster than standard MM. The order-up-to level that is set by any of the three methods is driven by the overestimation of λ (leading to overshooting the order level and fill rate) and the underestimation of μ c (leading to an undershoot).
For λ = 1 / 16 , μ c = 2 , and L = 2 the underestimation of μ c dominates for all three methods, leading to fill rates that converge from below. Although the overestimation of λ is more severe, this parameter affects only the inventory level distribution (via the lead time demand), whereas μ c affects also the distribution of a single customer's demand size, which both appear in the fill rate expression. Whereas the proposed method and ML are already within 0.1% from the 95% fill rate target for a sample of 50 observations, standard MM only attains 93.8% even for n = 200 .
When λ = 1 / 4 (see Fig. 3 b), the overestimation of λ and underestimation of μ c by the proposed method and ML almost cancel each other out at the order level calculation, leading to a realized fill rate very close to its target over the entire range of sample sizes. In Fig. 3 c ( λ = 1 ) all methods lead to fill rate overshoots, indicating that the overestimation of λ is dominant. ML performs only marginally better than the proposed method and the differences between the three approaches are smaller than in the previous scenarios. In the last scenario ( Fig. 3 d), λ is reset to 1/4, but the lead time is quadrupled to 8 periods. Comparing this scenario with that where λ = 1 / 4 and L = 2 , we observe that standard MM actually performs better under a longer lead time. By extending the lead time, more weight is given to the arrival rate, so that the underestimation of μ c is counteracted.
Summarizing all results, we find that the standard MM estimators show slow convergence of achieved fill rates, especially for highly intermittent demand patterns. ML and the proposed method perform closer to their targets. For low arrival rates, the proposed method, ML, and standard MM all converge to the optimal fill rate from below, whereas for higher arrival rates this pattern shifts to a convergence from above. The analysis in this section shows thatgiven that demands are stored per period -the best period length for optimal inventory performance is a non-trivial function of the arrival rate, demand size mean, and lead time.

Robustness analysis
In real forecasting and inventory control applications one should not only estimate demand parameters, but also select the distributional form. In this section we analyze the robustness of the different methods to misspecification of the demand size distribution. We first analytically demonstrate the robustness of the proposed method and then numerically compare its performance to that of standard MM and ML. Under the demand model specified in Section 2.1 , the proposed estimates for the arrival rate and the demand size mean are ˆ λ = − ln (n 0 /n ) and ˆ μ c = ˆ c is the variance of the compounding distribution [40] . Therefore, . Substituting μ c = 1 /λ and σ 2 c = 1 /λ 2 for the exponential distribution, the right hand side simplifies to λ, but if the compounding distribution is not exponential, then this does not hold true generally. Similarly, we can derive that For the exponential compounding distribution, the right hand side simplifies to 1 /λ ≡ μ c , but again this does not hold true generally. We remark that for other compounding distributions, the asymptotic properties of any estimator can be derived in a similar way. For the asymptotic analysis of the ML estimator we refer to Öztürk [30] . The robustness of the proposed estimator extends to compounding distributions with multiple parameters. The rationale is that it uses P central moments and the probability of zero demand for estimating P parameters of the compounding distribution and the arrival rate. The P th central moment of the full compound Poisson distribution depends only on the first P central moments of the compounding distribution [40] . This is not the case for the standard MM estimator, which always uses P + 1 moments for P parameters of the compounding distribution and the arrival rate. This makes the MM estimator sensitive to the (next moment of the) compounding distribution and therefore generally inconsistent under distribution misspecification.
We demonstrate this for the case where the compounding distribution has two parameters, which we take to be the mean μ c and the variance σ 2 c . The proposed estimator is now the solution to the following three equations ˆ λ = − ln (n 0 /n ) The first two equations directly lead to the same estimators for λ and μ c as for a one-parameter compounding distribution. Substituting these into the last equation gives This closed form again does not depend on the specific compounding distribution. Since this equation is directly obtained from the general formula of the second moment of a compound Poisson distribution, and as ˆ λ → λ and ˆ μ c → μ c , the continuous mapping theorem implies that ˆ σ 2 c → σ 2 c , irrespective of the compounding distribution. The standard MM estimator in this case replaces the first equation by the third central moment equation, which depends on the first three central moments of the compounding distribution [40] . Therefore, it uses the third moment of the compounding distribution to estimate its second moment. If more pa-rameters (or moments of the compounding distribution) are to be estimated, then the recursive solution procedure of the proposed estimator (and therefore also its robustness) readily extends.
We continue to numerically analyze the achieved accuracy of standard MM, ML, and the proposed estimator under misspecification of the compounding distribution. In line with the remainder of the paper, we assume that the demand sizes are exponentially distributed. However, we use two different distributions to generate demand sizes, namely the (continuous) uniform distribution with bounds 0 and 2 μ c , and the (discrete) logarithmic distribution (which was also used by Axsäter [2] ) with mean μ c . The proposed estimators for the arrival rate and the demand size mean have the same closed form for any compounding distribution and converge to the true values by construction. The other estimators, however, do show asymptotic biases under misspecification. Fig. 4 shows the results for the uniform distribution.
The standard MM estimator for λ shows an asymptotic upward bias of 52% in all scenarios, whereas for μ c it shows an asymptotic downward bias of 33% in all scenarios. This method therefore does not yield reliable estimates at all when the demand distribution is misspecified. ML is more robust than standard MM, but its asymptotic bias is still considerable for larger values of λ. For λ = 1 / 16 the error is negligible. For λ = 1 / 4 we record for n = 200 a 0.7% upward bias in λ and a 0.3% downward bias in μ c . For λ = 1 this error has increased to +5% for λ and -4% for μ c . The proposed method is, in line with our expectation, insensitive to the misspecified demand distribution and shows an unabated accuracy in comparison to the analysis of Section 4.2 . Also in the distribution   5 shows the results for the logarithmic distribution. As the logarithmic distribution is discrete, it can happen (especially for small sample sizes) that all observations have the same value and the variance is zero. The standard MM estimator (for the exponential compounding distribution) is then undefined, which is the reason why in Fig. 5 e and f the bias of standard MM is not drawn for n ≤ 7 . This issue does not exist for ML and the proposed estimator. The plots are qualitatively similar to those of Fig. 4 . As the logarithmic distribution's shape more closely resembles that of the exponential distribution, the asymptotic biases of ML and standard MM are smaller in this case. Nevertheless, standard MM still shows a large bias in all scenarios, and for λ = 1 we also see a clear asymptotic bias in ML. Only the proposed estimator is consistent (and has the smallest bias for all sample sizes) in all set-ups.

Case study
In this section we compare the different estimators for compound Poisson demand parameters in practice, using a data set of demands for 496 different commercial airline parts. This data set contains demand observations at the individual customer order level, which we aggregate to 125 weekly totals. We restrict our analysis to the 404 parts that have had at least one demand in the first 20 weeks. In Table 4 we present some descriptive statistics. The estimated arrival rates and demand size means are calculated based on the individually stored demands of the full history per item. Almost all arrival rates are below 1, whereas the demand sizes vary wildly.
Based on the weekly aggregated demand observations, we set up a comparative analysis for the arrival rates and demand size means according to standard MM, ML, and the proposed method. As a benchmark, we furthermore include the 'theoretically best' estimates based on the individual orders, which, of course, cannot be exploited using aggregated data. Lengu et al. [12] established that a compound Poisson process with geometric compounding distribution has a good fit to this data. We adopt this distribution (which is further studied in Appendix B ), but in line with the rest of this paper, we also test the exponential compounding distribution. Starting with the first 20 weeks and every time increasing the time horizon by one week, we estimate for every horizon the arrival rate and demand size mean for all items. We report the Mean Absolute Percentage Error (MAPE) between these estimates and the 'theoretically best' estimator over the entire sample (for every time horizon, as an average across all items). The proposed method and the 'theoretically best' method do not depend on the specific compounding distribution. Furthermore, the ML results do not vary visibly between the exponential and geometric compounding distribution. Therefore, only for standard MM do we report the results under both distributions separately. The results are depicted in Fig. 6 . By construction, the MAPE of the individual demands estimator converges to 0 and the evolution of this estimator throughout the sample benchmarks the other methods. The standard MM estimator clearly performs worse when the exponential distribution is assumed. However, in both cases it performs considerably worse than the other methods. For the arrival rate, the proposed estimator slightly outperforms ML by 0.2% for the full sample. For smaller samples, the differences between the estimators are smaller, but their performance ranking does not change. For the demand size mean, standard MM even deteriorates when the sample size passes 90 weeks. The theoretical results from the preceding sections are confirmed also in this practical case study: standard MM performs by far worst and is most sensitive to the exact specification of the compounding distribution, and the proposed estimator marginally outperforms ML. It is furthermore established that the accuracy loss due to storing only weekly data increases with the sample size. If all 125 weeks are used, then the proposed estimator loses approximately 10% in accuracy compared to using individual transaction data. Standard MM, contrarily, loses approximately 35% or 40%, depending on the compounding distribution.

Conclusions
The inventory control literature provides hardly any guidance on obtaining demand parameters from period demand data, although companies typically do store historical demand observations per period and software packages also use period demand as input for their forecasts. Period demand forecasts cannot be directly transformed into compound Poisson parameter estimates, even if the forecasting procedure yields separate estimates for the time between periods with positive demands and the average pe-riod demand size. Misusing these to obtain demand parameters leads -also asymptotically -to severely biased estimates, overshot fill rates, and thus excessive inventories and associated costs. It is important to observe that this flaw occurs if any period forecasting procedure is used to estimate individual customer demand parameters. We have addressed this lack of connection between the forecasting and inventory control literature, and performed a comparative study of different compound Poisson parameter estimators. We proposed an alternative MM estimator that uses a nonstandard choice of moments, reflecting the intermittent nature of the demand. The fraction of periods without demand is used to estimate the arrival rate, and the estimated arrival rate is combined with further central moments of the (full) compound Poisson distribution to estimate the parameters of the compounding distribution (the distribution of the individual demand sizes).
Our results confirmed that the standard MM estimator (the default option in inventory control) has severe biases for a wide range of intermittent demand patterns and sample sizes, and has a larger variance than the other studied methods. The proposed method's deviation from the true value of both parameters is smaller over the entire range of intermittent demand patterns and sample sizes considered. In several scenarios (particularly also in our case study) the proposed estimator has a lower bias than ML, whereas in the other scenarios it performs almost indistinguishably. A key advantage of the proposed method is that its closed form does not depend on the specific demand size distribution, which implies that it is completely robust to misspecification of that distribution. This property (which neither standard MM, nor ML possesses) is particularly relevant in practice. The closed-form nature of the estimator is useful in large-scale applications where many items are controlled simultaneously and computation speed is an issue. Both practitioners and applied researchers should be aware of 1) the accuracy loss due to storing demand only periodically, and 2) the consequences of their estimator choice on accuracy and inventory performance. It is furthermore important to realize that this parameter estimation issue exists also for other inventory models and service measures, as well as outside inventory control. Accuracy results carry over directly, but the final effect is of course specific to the decision model.
An important empirical research avenue is to study the effect of bucketing the demand observations into different period lengths. Our study showed that the relationship between the intermittency of the demand pattern and the inventory performance of the different estimators is non-trivial. A further research opportunity lies in the fact that in practice it is sometimes questionable if and for how long the demand process remains stationary. A non-stationary compound Poisson process -which does not have a fixed arrival rate, but a time-varying arrival intensity function -can model this situation. If it is assumed that the intensity function is piece-wise linear, then the process reduces to a stationary compound Poisson process on the corresponding subsets of the data and our method can still be used. For more general intensity functions, however, existing literature is solely based on ML, yielding an opportunity for future research. Finally, an overall connection is lacking in the broader sense between the forecasting and inventory control literature. Forecast accuracy is typically assessed via symmetric loss measures, whereas from an inventory control perspective over-and underestimation carry asymmetric penalties. A study comparing different forecasting methods from this perspective may lead to new insights. For the ML estimators, we use that the probability mass function of the Poisson-geometric distribution is [32] f PG (x ) = exp (−λ) and follow the procedure of Section 3.2 . The proposed method does not depend on the compounding distribution and stays unaltered: ˆ λ = − ln (n 0 /n ) and ˆ μ c = −ˆ μ 1 ln (n 0 /n ) . We perform the same numerical experiment as outlined in Section 4.2 , now using the geometric compounding distribution both for generating the individual demands and calculating the estimates. We discuss only the setting with λ = 1 / 4 and μ c = 2 and omit the variance and fill rate analysis, as all results are very similar to those under the exponential distribution. Fig. 7 shows the results. Standard MM is now more accurate than under the exponential distribution, but still leaves a considerable gap with the other methods.