Robust compound Poisson parameter estimation for inventory control

University of Groningen

Robust compound Poisson parameter estimation for inventory control

Prak, Dennis; Teunter, Ruud; Babai, Mohamed Zied; Boylan, John E.; Syntetos, Aris

Published in:

Omega

DOI:

10.1016/j.omega.2021.102481

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from

it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2021

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Prak, D., Teunter, R., Babai, M. Z., Boylan, J. E., & Syntetos, A. (2021). Robust compound Poisson

parameter estimation for inventory control. Omega, 104, [102481].

https://doi.org/10.1016/j.omega.2021.102481

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

ContentslistsavailableatScienceDirect

Omega

journalhomepage:www.elsevier.com/locate/omega

Robust

compound

Poisson

parameter

estimation

for

inventory

control

R

Dennis

Prak

a,b,∗

_,

_Ruud

_Teunter

a

_,

_{Mohamed Zied}

_Babai

c

_,

_{John E.}

_Boylan

d

_,

_Aris

_Syntetos

e a Department of Operations, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

b Department Industrial Engineering and Business Information Systems, University of Twente, PO Box 217, 7500 AE Enschede, The Netherlands c Kedge Business School, 680 Cours de la Libération, 33405 Talence Cedex, France

d Department of Management Science, Lancaster University, Lancaster LA1 4YX, United Kingdom

e PARC Institute of Manufacturing, Logistics and Inventory, Cardiff Business School, Cardiff University, Cardiff CF10 3EU, United Kingdom

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 12 June 2020 Accepted 4 May 2021 Available online 6 May 2021

Keywords:

Inventory control Forecasting Probability

Compound Poisson demand

a

b

s

t

r

a

c

t

Mostcompaniesstoredemanddataperiodicallyandmakeperiodicdemandforecasts,whereasmany de-mand processesininventorycontrolneed parameterestimates atthe individualcustomerlevel. Guid-anceonestimatingtheparametersofacontinuous-timedemandprocess fromperioddemanddatais lacking,inparticularforthepopularandwell-studiedcompoundPoissonclassofdemand.Whereasthe statisticsliteraturetypically focusesonasymptoticproperties,parametersforinventorycontrolhaveto beestimatedbasedonalimitednumberofperiodichistoricaldemandobservations.Weshowthatthe standardMethod-of-Moments(MM)estimator– thedefaultchoiceinappliedinventorycontrolresearch – isseverelybiasedforfinitesamples.TheMaximumLikelihood(ML)estimator– whichneedstobe ob-tainedbyanumericalsearch– performsbetter,butbothestimatorslackrobustnesstomisspecificationof thedemandsizedistribution.Weproposeanintuitive,consistent,closed-formMMalternativethat dom-inatesstandardMM andML intermsofestimationaccuracyand on-targetinventoryperformance.Its closedformdoesnotdependonthespecificdemandsizedistribution,makingitrobustandeasily appli-cableinlarge-scaleapplicationswithmanyitems.Inacasestudy,wefindthattheaccuracylossdueto storingdemandperiodicallyisfourtimesashighunderstandardMMasundertheproposedestimator.

© 2021TheAuthor(s).PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

ThecompoundPoissondemandprocessisoneofthemost gen-eral and widely studied demand processes for inventory control. It has a natural business interpretation with customers arriving incontinuoustime,andtheirdemandsizesfollowinganygeneral compounding distribution, and is therefore a standard choice in textbooks (e.g.[1–3]) andcommercialsoftware(e.g.Forecast Pro, Oracle NetSuite, SAP APO, SAS, and Slimstock), especially for in-termittent demand patterns. After Galliher et al. [4] studied the specialcaseofPoisson-geometricdemand,FeeneyandSherbrooke

[5] and Chen etal.[6] derived order-up-topolicies undergeneral compoundPoissondemand.ArchibaldandSilver[7] extendedthis work to general

(

s,S

)

policies.Order-up-to policies canbe deter-R Area - Supply Chain Management This manuscript was processed by Associate

Editor Keskin.

∗_{Corresponding author. Present address: Department Industrial Engineering and}

Business Information Systems, University of Twente, PO Box 217, 7500 AE Enschede, The Netherlands

E-mail addresses: d.r.j.prak@utwente.nl (D. Prak), r.h.teunter@rug.nl (R. Teunter),

mohamed-zied.babai@kedgebs.com (M.Z. Babai), j.boylan@lancaster.ac.uk (J.E. Boy- lan), SyntetosA@cardiff.ac.uk (A. Syntetos).

mined efficientlyundercompound Poissondemand,asshownby Sherbrooke [8] and Graves [9], who provided approximate pro-cedures,and Axsäter [10], who presented an exact method. Sev-eralempiricalstudiesconcludedthat,especiallyifproductsarenot fast-moving, compound Poissondistributions provide the best fit (e.g.[11–13]).Therefore, they remain astandard choice in inven-torycontrolresearch(e.g.[14–16]).

In orderto exactly characterizeinventory policies, a complete specificationofthecompoundPoissondemandprocessisrequired. This is most obvious for continuous-review inventory systems. Here,thearrivalofacustomer immediatelytriggersan order,and the probability distribution of individual customer demand sizes determines the expectedreorder level undershoot. However, also underperiodicreviewallparametersofthecompoundPoisson de-mandprocessareneededforexactcalculations[17].Ascompound Poisson demand processes operate in continuous time, it would be natural to estimate their parameters fromindividual transac-tiondata.The averagenumberofcustomer arrivalsin aperiodis anunbiasedestimatorforthearrivalrateandtheparameter(s)of thecompounding distribution canbe estimateddirectly fromthe demand sizes. However, it is common in practice (and a typical assumption inapplied inventorycontrol literature [12,13,18])that

https://doi.org/10.1016/j.omega.2021.102481

only total demand per period is stored and used as a basis for forecasting.Infact,temporaldataaggregationtoimproveaccuracy is frequentlysuggested(see e.g.[19–21]). Asa result, commercial forecasting software is mainly based on period demand. Storing demandperperiodratherthanperindividualtransactiondoeslead to anaccuracyloss.As weshow inthispaper,boththechoice of theparameterestimatorandtheperiodlengthtowhichdemands are aggregated have a large influence on the magnitude of the error.

Separate forecastingofdemandsizesanddemandinter-arrival times may be achieved by so-called Size-Interval methods. The most well-known one – twice applying exponential smoothing– wasfirstproposedbyCroston[22].Asurveyofsoftwarepackages’ forecasting capabilities [23] confirms the popularity of Croston’s method. Continuations of this research stream are reviewed by Syntetosetal.[24].However,periodforecastingmethodsestimate themeandemandsizeinaperiodandthemeanintervalbetween periods withdemand occurrence, rather than the mean demand sizeofanindividualcustomerandmeanintervalbetween individ-ual customer arrivals. Therefore, they cannot be directly used in inventorycontrol,especiallynotifdemandisintermittent[25,26]. Ignoring thismismatch leadseven asymptotically to significantly biasedestimatesanddramaticallyunderachievedservicelevels,as wewillshowinSection2.

Very few authorsin theinventory control literature discussed howtheunknown parametersofcompoundPoissondemand pro-cesses (or ofdemandprocesses ingeneral)can be obtainedfrom period demand data. The authors who did,suggested the use of standard Method-of-Moments (MM) estimators [2,27], which are thereforeatpresentthestandardchoiceinappliedwork[13,18].In the wider (statistics) literature, the discussion of parameter esti-mationcentersaround(variantsof)MMandMaximumLikelihood (ML). Examples are Anscombe [28] and Sichel [29] (logarithmic compoundingdistribution),Öztürk[30] (exponentialcompounding distribution), Withers and Nadarajah [31] (gamma compounding distribution), and Balakrishnanet al. [32] (geometric compound-ing distribution). MM and ML variants are also givenfor certain compoundPoissondistributionsintextbookssuchasJohnsonetal.

[33],withoutaperformancecomparison.Ifauthorscompare prop-ertiesofdifferentestimators,thenthediscussionisalmostalways limitedtolarge-sampleproperties.Inthispaperwefocuson finite-sampleperformanceinstead,whichisimportantinpractice. Prod-uctlife-cyclesareshortening,andevenifhistoricdemanddataare availableformultipleyears,thennotallofthemmaybeindicative offuturedemand.

OurmaincontributionisthepresentationofanalternativeMM estimator that retains Croston’s core idea of separately estimat-ing customer arrival rates and demand size means. It bases the customer arrival rateon the fractionofperiods without demand, andthenusescentralmomentstoestimate theparameters ofthe (individual) demand sizedistribution.This estimatorhasan intu-itive,closed-formsolution,independentofthespecific compound-ingdistribution.Itisfurthermorerobusttomisspecificationofthe compounding distribution,whichisanimportantpractical advan-tage. Theideaofestimatingdistributionparameters basedonthe ‘zerofrequency’hassporadicallybeenapplied[12,28,34–36]. How-ever, i) ageneral methodto estimate anycompound Poisson de-mandparametersbasedonthislogichasneverbeengiven,andii) astudyofsuchamethod’sperformanceislacking.

Focusing on the exponential compounding distribution in the main bodyofthepaperanddiscussingthegeometric compound-ing distributioninAppendixB,wecomparethefinite-samplebias and variance of standard MM, ML, and the proposed alternative MM estimator.Thereafter,we comparetheachievedfillratesina continuous review order-up-to system. We furthermore study all estimators undermisspecificationofthe demandsizedistribution

andfindthatbothstandardMMandMLlackrobustness, whereas the proposed method is unaffected by the specific demand size distribution.A casestudyvalidates theresults andquantifiesthe 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 servicemeasure, and discusses the misapplication of period Size-Intervalmethods.Section3 introducestheconsistentstandard MMandMLestimatorandtheproposed estimator.Section 4 pro-vides a comparison of the consistent estimators’ biases, whereas

Section 5 discusses their variances. Section 6 presents a com-parison of achieved fill rates in an order-up-to inventory model, and Section 7 gives robustness results. A comparison of the estimators in a realcase studyand a quantification ofthe accu-racy loss due to storing aggregated demands per period rather than individual transactions is givenin Section 8.Section 9 con-cludes. Appendix A contains derivations related to Section 4, andAppendix B presents results forthe geometric compounding distribution.

2. Demandmodel,inventorymodel,andsize-intervalmethods In this section we specify the general compound Poisson de-mandmodelandthecompoundingdistributionsonwhichwebase ourresults.Thereafterwediscusstheorder-up-toinventorymodel andthefillratecalculation.Finally,wediscusstheconsequencesof misusingperiod Size-Intervalestimates ascompoundPoisson de-mandparameters.

2.1. Demandmodel

Weconsideramodelwheredemandfollowsacompound Pois-sonprocesswithsomegeneral(discreteorcontinuous) compound-ing distribution. We have a sample of n period demand realiza-tionsx1,...,xn.Thatis,fort=1,2,...,n,weobserverealizationsof therandom variableXt=Ki=1Di,whereK follows aPoisson dis-tributionwithmean

λ

and,independently,theindividualcustomer demands Di are i.i.d.withmean

μ

c. Theindividual customer de-mandsare not observed. Werequirethat P

(

Di≤ 0

)

=0asa cus-tomercannaturallyonlyhavepositivedemand.Thesetx1,...,xnof perioddemandrealizations isusedtoestimate theparameters of thecompoundPoissonprocess.

In the remainder of the paper we focus on the exponential compoundingdistributionwhichiscontinuousandthereforemore suitablefordemonstration,sinceanyfillratelevelcanbeachieved exactly.WerefertoAppendixB foradditionalresultsbasedonthe geometriccompoundingdistribution.Fortheexponential distribu-tion,the individual demand sizesDi have theprobability density function

fE

(

x

)

=

θ

exp

(−

θ

x

)

,

forx_{≥ 0}and

θ

_>0,withdemandsizemean

μ

c=1/θ.

Forthegeometric distribution,the individual demandsizesDi havethe(discrete)probabilitymassfunction

fG

(

x

)

=

(

1−

β

)

β

x−1,

forx₌1_,2_,...and0_<

β

_<1,withdemandsizemean

μ

c=1/

(

1−

β

)

.

2.2. Inventorymodel

We study a continuous review order-up-to inventory control modelwitha positivefixed lead time ofLperiods andfull back-ordering,underafillrateconstraint.Order-up-toinventorymodels

are relatively easy to optimize, popularin practice,and theoreti-callyoptimalforvariouscostandservicemodelswhenordercosts arenegligible[1,37].Asaservicetargetweusethefillrate,which istheexpectedfractionofademandthatcanbeserviced immedi-atelyfromon-handinventory.

We followthecomputationofthefillrateforcompound Pois-sondemandmodelsbyAxsäter[2],andgeneralizeittohandlealso continuouscompoundingdistributions.Specifically,denotingan in-dividual demandby Di andtheinventory level by IL, fora given order-up-tolevelS,thefillrateistheexpectedfractionofasingle orderthatcanbefulfilledfromon-handstock,

FR

(

S

)

=EIL



EDi[min

(

max

(

IL,0

)

,Di

)

]



μ

c .

ThedistributionofDiistheselectedcompoundingdistribution withmean

μ

c,whereasthedistributionofILfollowsfromthe ob-servation that the inventory level attime t₊L isequal to S mi-nustheleadtimedemand.Sincedemandperperiodiscompound Poissonwithrate

λ

anddemandsizemean

μ

c,thedemandduring the lead time follows a compound Poissondistribution withthe same compounding distribution and rate

λ

L. Denoting lead time demand by DL, we findthe relationshipP

(

IL≤ x

)

=P

(

DL≥ S− x

)

. So,FIL

(

x

)

=1− FD_L

(

S− x

)

,whereFIL denotes thedistribution func-tionoftheinventorylevelandFD_Ldenotesthedistributionfunction ofthecompoundPoissondistributionofthedemandprocess,with arrivalrate

λ

Landdemandsizemean

μ

c.Theobjectiveis

minimize S

subjectto FR

(

S

)

≥

γ

,

where

γ

isthetarget fillrate. Intheremainder ofthispaper,we use

γ

=0.95(i.e.95%).

2.3. Misapplicationofperiodsize-intervalmethods

Period Size-Interval methods yield estimates of the time be-tween two periods withpositive demand andthe size ofa non-zero period demand. We demonstrate here the consequences of misapplyingtheseascompoundPoissonparameterestimators.We firstdiscussCroston’smethodanda‘bias-corrected’versionofthis method. Thereafterwe discuss an intuitive method that usesthe average of all non-zero period demands and the average of all timesbetweenperiodswithpositivedemand,whichwerefertoas ‘UnweightedAveraging’.Alldiscussedmethodsonlyupdateaftera periodwithpositivedemand.

Croston’s methodestimatesinperiodnthe timebetweentwo 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 haveobserved n periods,the nthperiodhasapositivedemandxn,andthepreviouspositive de-mandwasqnperiodsago,

ˆ

sn=

α

1xn+

(

1−

α

1

)

sˆn−1,

ˆ

pn=

α

2qn+

(

1−

α

2

)

pˆn−1,

where0<

α

1<1and0<

α

2<1areconstants.

IfCroston’speriodestimatesarewronglytakenastheestimates oftheindividualdemandparameters,thensˆnistheestimate after n period observations forthe mean

μ

c ofthe compounding dis-tribution, whereas 1_/pˆn is the estimate for

λ

. As shownby Syn-tetos and Boylan [38], using 1/pˆn as an estimate for the arrival rate of a compound Poisson process leads to an inversion bias, asE[1_/pˆn]=1/E[pˆn].Thisinversionbiasdoesnotdiminishasthe sample size goes to infinity, as Croston’s method uses exponen-tial smoothingandthus doesnot useall datapoints, so thatthe CentralLimitTheoremdoesnothold.SyntetosandBoylan[39] ap-proximatelycorrectedforthisbiasusingasecond-orderTaylor ex-pansion, andderived theestimator

(

1−

α

2/2

)

/pˆn for

λ

.The

esti-matorfor

μ

cdoesnotsufferfromtheinversionbiasandremains unaltered.

Itisstraightforwardtoshowthat,withorwithouttheinversion biascorrection,Croston’s methodyieldsinconsistentcustomer de-mandparameterestimateswhen appliedtoperioddemand. Con-sidersˆn,whichisunaffectedbytheinversionbias.Insteadof esti-matingthesizeofanindividualdemand,itestimatesthesizeofa positiveperioddemand.Therefore,

lim n→∞E[sˆn]=E

(

Xn

|

Xn>0

)

=

μ

c _∞ k=1kexp(−λ)λ k k! 1− exp

(−

λ

)

= 1− exp

μ

c

λ

(−

λ

)

, whichisstrictlylargerthan

μ

c andsmallerthan

μ

c

λ

forall

λ

>0 and

μ

c>0.Thislimitapproaches

μ

cas

λ

→0,and

μ

c

λ

as

λ

→∞. Theestimator’slimit(andthusitsbias)isproportionalto

μ

c.

Theasymptoticbiasofthearrivalrateestimateismore cumber-sometoquantify.Insteadofestimatingthetimebetweentwo cus-tomer arrivals, pˆn estimates thenumber ofperiods betweentwo periods with positive demand. This number of periods (denoted by q) follows a geometric distribution with success probability 1− exp

(

−

λ

)

, mean 1/

(

1− exp

(

−

λ

))

, and variance exp

(

−

λ

)

/

(

1− exp

(

₋

λ

))

2_{.AlthoughE[pˆ}

n]→1/

(

1− exp

(

−

λ

))

asn→∞,we can onlyapproximatethelimitofE[1/pˆn].Wedosobyasecond-order Taylorexpansionandfind,usingthemeanandvarianceofthe ge-ometricdistributionandthefactthatpˆnisanexponential smooth-ingestimatorwithasymptoticvarianceVar(q

)

α

2/

(

2−

α

2

)

,that

lim n→∞E[1/pˆn]≈ 1 E[q]+

α

2 2−

α

2 Var

(

q

)

E[q]3 =



1+

α

2 2−

α

2 exp

(−

λ

)



(

1− exp

(−

λ

))

.

As

α

2→0, this limit approaches 1− exp

(

−

λ

)

, which is the

re-sultwithoutinversion bias.Inthat scenario,

λ

isstrictly underes-timatedforall

λ

>0.Theapproximatelimitoftheestimator con-vergesto0as

λ

→0,andto1as

λ

→∞.Obviously,thebias cor-rectionproposedbySyntetosandBoylan[39] scalestheestimated arrivalrate,andthusitsexpectation,by

(

1₋

α

2/2

)

.Inconclusion,

Croston’smethodgenerallyunderestimatesthearrivalrate,butthis issomewhatcounteractedbyitsinversionbias,leadingtoan over-estimationespeciallyforlow valuesof

λ

and/orhighvaluesof

α

2.

However,(approximately)correctingforthat inversionbias inten-sifiestheunderestimation.

WhileCroston’smethodanditsvariantsgive greaterweightto themostrecentdemandsizesandintervals,thisisnota require-mentforSize-Interval methods.An alternativeapproachis touse equalweights, butretain theseparate estimation ofdemand size meansandintervals. TheUnweightedAveragingmethoddoesthis bytakingasthedemandsizemeantheaverageofallpositive pe-riod demand sizes, andasarrival rate theinverse ofthe average number of periods between two periods with positive demands. Giventhatwehaveobservednperiods,ofwhichnphadapositive demand,then−thperiodhasa positivedemandxn,andwehave observedthepreviouspositivedemandqnperiodsago,Unweighted Averagingsets ˆ sn= 1 np xn+ np− 1 np ˆ sn−1, ˆ pn= 1 np qn+ np− 1 np ˆ pn−1.

Itisstraightforwardtoobservethat,justaswithbothCroston vari-ants,

lim

n→∞E[sˆn]=E

(

Xn

|

Xn>0

)

=

μ

c

λ

1− exp

(−

λ

)

.

Unlike Croston’s method, the Unweighted Averaging method is asymptotically unbiased, since it uses all positive demands and

Table 1

Asymptotic Biases and Achieved Fill Rates when Using Period Size-Interval Methods for Compound Poisson Parameter Estimation (Target Fill Rate 95%, Lead Time L = 2 ).

Achieved Fill Rates (Compounding Distr.) True Parameters Asymptotic Bias (Exponential) (Geometric)

λ μc α2 λCROST λSBA λUA μc Standard SBA UA Standard SBA UA

1/16 2 0.1 + 1.7% −3.4% −3.1% + 3.2% 95.5% 95.4% 95.4% 97.2% 97.2% 97.2% 1/16 2 0.3 + 13.1% −3.9% −3.1% + 3.2% 95.6% 95.4% 95.4% 97.2% 97.2% 97.2% 1/16 2 0.5 + 27.3% −4.5% −3.1% + 3.2% 95.8% 95.4% 95.4% 97.2% 97.2% 97.2% 1/16 5 0.1 + 1.7% −3.4% −3.1% + 3.2% 95.5% 95.4% 95.4% 95.7% 95.7% 95.7% 1/16 5 0.3 + 13.1% −3.9% −3.1% + 3.2% 95.6% 95.4% 95.4% 96.5% 95.7% 95.7% 1/16 5 0.5 + 27.3% −4.5% −3.1% + 3.2% 95.8% 95.4% 95.4% 96.5% 95.7% 95.7% 1/4 2 0.1 −7.9% −12.5% −11.5% + 13.0% 96.5% 96.3% 96.3% 97.0% 97.0% 97.0% 1/4 2 0.3 + 0.6% −14.5% −11.5% + 13.0% 96.8% 96.2% 96.3% 98.2% 97.0% 97.0% 1/4 2 0.5 + 11.4% −16.4% −11.5% + 13.0% 97.1% 96.1% 96.3% 98.2% 97.0% 97.0% 1/4 5 0.1 −7.9% −12.5% −11.5% + 13.0% 96.5% 96.3% 96.3% 97.0% 96.4% 97.0% 1/4 5 0.3 + 0.6% −14.5% −11.5% + 13.0% 96.8% 96.2% 96.3% 97.0% 96.4% 97.0% 1/4 5 0.5 + 11.4% −16.4% −11.5% + 13.0% 97.1% 96.1% 96.3% 97.4% 96.4% 97.0% 1 2 0.1 −35.6% −38.8% −36.8% + 58.2% 98.6% 98.4% 98.5% 99.3% 99.0% 99.3% 1 2 0.3 −32.7% −42.8% −36.8% + 58.2% 98.7% 98.1% 98.5% 99.3% 99.0% 99.3% 1 2 0.5 −29.0% −46.8% −36.8% + 58.2% 98.9% 97.8% 98.5% 99.3% 98.5% 99.3% 1 5 0.1 −35.6% −38.8% −36.8% + 58.2% 98.6% 98.4% 98.5% 98.8% 98.6% 98.8% 1 5 0.3 −32.7% −42.8% −36.8% + 58.2% 98.7% 98.1% 98.5% 98.9% 98.4% 98.8% 1 5 0.5 −29.0% −46.8% −36.8% + 58.2% 98.9% 97.8% 98.5% 99.1% 98.2% 98.8%

inter-arrival timesin the dataset ratherthan applying diminish-ingweightsviaexponentialsmoothing.Therefore,theCentralLimit Theoremholdsandpˆn convergesinprobabilitytoE[pˆn],sothat

lim

n→∞E[

λ

ˆn]=nlim→∞E[1/pˆn]=1− exp

(−

λ

)

.

The UnweightedAveragingestimatorfor

λ

convergesto 0as

λ

→ 0, and to 1as

λ

→∞.This reveals the flaw in using period de-mandestimatorsasindividualdemand parameters:nodistinction canbemadebetween1orseveraldemandsinaperiod,leadingto an underestimatedarrivalrate.Allindividualdemandsinaperiod areaddedandconsideredasasinglecustomerdemand,leadingto anoverestimateddemandsizemean.

Table 1 quantifiesthe (approximate)asymptotic biasresulting fromusingstandardCroston’smethod(‘CROST’),thecorrected ver-sionbySyntetosandBoylan[39] (Syntetos-BoylanApproximation, ‘SBA’),ortheUnweightedAveragingmethod(‘UA’),forseveral pa-rameter choices. Furthermore, we report the fillrate that would be achievedforboththeexponentialandgeometriccompounding distribution,underthetruedemandparameters,whenatargetfill rateof95%ischosen.Weusealeadtimeof2periods.

None oftheresultsdepend on

μ

c,exceptforthe achievedfill ratesunderthegeometriccompoundingdistribution.Thereasonis thatfordiscretedemanddistributions,fillratescannotbeachieved exactly. For low arrival rates and/or high smoothing parameter values, Croston’s inversion bias results in an overall overestima-tionof

λ

.TheSyntetos-BoylanApproximationachievesan asymp-toticunderestimation forallscenarios.The UnweightedAveraging methoduniformlyunderestimates

λ

.Thedemandsizemean

μ

cis alwaysoverestimatedforallmethods,andall(absolute)percentage asymptoticbiasesareincreasingin

λ

.

Interestingly,all achievedfillratesareasymptotically toohigh, irrespective ofwhether

λ

is over- or underestimated. Thisis be-cause the achievedfill ratedependsmore strongly on

μ

c, aswe willfurtherdiscussinSection6.Theachievedfillratesare increas-ing in

λ

, andfor

λ

=1 they are over 97.8%for all compounding distributions. We conclude that misapplying period Size-Interval methodstoestimate compoundPoissondemandparametersleads to severely overestimated demand size means, significantly over-shotfillratesandcorrespondinglyexcessiveinventorycosts.

3. Consistentestimators

Having discussed the misapplication of period forecasting methodsforobtainingcompoundPoissondemandparameters,we movetoconsistentestimationmethods.Werevisit thecompound Poissondemandmodelintroducedintheprevioussectionand dis-cussthe standardMMestimator(appliedintheinventory control literature),theMLestimator,andtheproposedalternativeMM es-timator - first for general compound Poisson demand processes andthenforthespecificcaseofanexponentialcompounding dis-tribution.

WefirstgivesomepropertiesofthePoisson-exponential distri-butionthatarehelpfulforouranalysis[30].Itsprobabilitydensity functionisgivenby fPE

(

x

)

=exp

(−

λ

− x/

μ

c

)

∞  i=1



(

λ

/

μ

c

)

ixi−1 i!

(

i− 1)!

,

for

λ

>0,

μ

c>0,andx>0.Thisdistributionhasa positive prob-ability massof exp

(

−

λ

)

atx=0.The first four central moments are

μ

1=

λμ

c,

μ

2=2

λμ

2c,

μ

3=6

λμ

3c,and

μ

4=

(

24

λ

+12

λ

2

)

μ

4c. Weassumethat asamplex1,...,xn ofperioddemands isavailable anddenotetheithcentralsamplemomentoftheperioddemand observationsby

μ

ˆi.Inparticular,

μ

ˆ1 denotesthesamplemeanand

ˆ

μ

2thesamplevariance.

3.1. Thestandardmethod-of-momentsestimator

StandardMM estimators equate central momentsof the com-poundPoissondemanddistributiontotheircorrespondingsample equivalents.Thatis,ifthecompoundingdistributionhasP param-eters,wesolve

μ

1 =

μ

ˆ1, . . .

μ

P+1 =

μ

ˆP+1,

for the unknown arrival rate and the P parameters of the com-poundingdistribution.

An exponential compounding distribution has one parameter, whichwe take to be its mean.So,there are intotal two param-eterstobeestimated. ThestandardMM estimatorsforthearrival

rateandthedemandsizemeanarethen ˆ

λ

=2

μ

ˆ21 ˆ

μ

2 and

μ

ˆc= ˆ

μ

2 2

μ

ˆ1.

MMestimatorsareguaranteedtobeconsistent,butmaystillhave severefinite-samplebiasesandvariances, asournumericalresults will reveal. Furthermore, the second order sample moment (the variance of the compound Poisson distribution) depends on the exact shape ofthe compoundingdistribution andnot onlyon its mean.ThisimpliesthatstandardMMisnotrobusttomisspecified demandsizedistributions.

3.2. Themaximumlikelihoodestimator

The MLestimatorofacompoundPoissondistributionisfound by maximizingthelog-likelihoodfunctionwithrespecttothe un-knownparameters.Inmostcases,themaximizing arguments can-notbefoundinclosedform,andanumericalsearchprocedurehas tobeinvoked.Thelog-likelihoodfunctionisdefinedas

L

(

θ|

x1,...,xn

)

= n  t=1

ln

(

fCP

(

xt

|

θ

))

,

where fCP denotes the probability density function of the corre-sponding (full) compound Poissondistribution and

θ

denotes its vectorofparameters,includingthearrivalrateofthePoisson pro-cessandtheparametersofthecompoundingdistribution.

For the exponential distribution, we have fCP = fPE and

θ

=

(

λ,

μ

c

)

. One can maximize L

(

λ,

μ

c

|

x1,...,xn

)

to find estimates

λ

ˆ and

μ

ˆc,using anon-linear maximizationalgorithm. AlthoughML iscomputationally moredemandingthanthe othermethods,itis not prohibitively expensive,asitinvolvesonly atwo-dimensional search.Moreover,itisguaranteedtoprovideaconsistentestimator thatasymptoticallyhasthelowestvariance.However,MLoffersno guaranteedoptimalperformanceforfinitesamples.

3.3. Theproposedmethod

We propose an alternative MM estimator andchoose its mo-ments such that the intuition of the Size-Interval forecasting lit-erature– separateestimationofthecustomerarrivalrateandthe demandsizemean– isemployed.Alsothisestimatorisguaranteed tobeconsistent,asitisanMMestimator.Wefirstdiscuss estima-tion ofthe arrival rate

λ

. Denote by N0 the (stochastic) number

ofperiods inthesampleofsizeninwhich nodemandoccurred. Thecorrespondingrealizationofthatnumberofperiodsisdenoted by n0.The occurrenceofnodemand inaperiod isBernoulli

dis-tributed with probability exp

(

₋

λ

)

, and therefore the number of periods withoutdemandout ofntotalperiodsfollowsabinomial distributionwithparametersnandexp

(

−

λ

)

.Therefore

E[N0]=nexp

(−

λ

)

.

Weequatethissamplemomenttoitsexpectationtofindthe esti-matorforthearrivalrate:

n exp

(−ˆ

λ

)

=n0,

leadingto ˆ

λ

=− ln

(

n0/n

)

.

For each of the parameters of the compounding distribution, we equate(aswithstandardMM)centralmoments ofthecompound Poissondistributiontotheir sampleequivalents.However, the ex-plicit use of the fractionof periods without demand to estimate thearrivalrateimpliesthatthehighestorderofmomentsusedby

this method is always one lower than for the standard MM ap-proach. Thatis, if the compounding distribution hasP unknown parameters,wesolve ˆ

λ

=− ln

(

n0/n

)

,

μ

1 =

μ

ˆ1, . . .

μ

P =

μ

ˆP.

Inthecaseofan exponential compoundingdistribution(more generally,inthecaseofanycompoundingdistributionwitha sin-gleparameter),thesetofequationsreducesto

ˆ

λ

=− ln

(

n0/n

)

,

μ

1 =

μ

ˆ1.

Using the fact that

μ

1=

λμ

c, we solve for the demand size mean:

ˆ

λ

μ

ˆc=

μ

ˆ1

where

μ

ˆ1 isthesampleaverageoftotaldemand perperiod.This

leadsto ˆ

μ

c=−

μ

ˆ1 ln

(

n0/n

)

.

Tocalculatetheseestimates,oneonlyneedstoobtainthe frac-tionofperiodswithoutdemandandtheoverallmeandemand dur-ing all periods.Furthermore,whereas standard MMandML have tobe derived separately forevery compounding distribution,this methodgives estimatesthat do notdepend on thespecific com-pounding distribution. Irrespective of the compounding distribu-tion,

λ

ˆ and

μ

ˆc have exactly the same closed-form solution. In

Section 7 we elaborate on this property, which extends to com-pounding distributions with an arbitrary number of parameters. Furthermore,we showthatthispropertyimpliesrobustness: irre-spectiveofthecompoundingdistribution,theproposedestimators convergetothecorrespondingtrueparameters.

4. Biascomparison

Inthis section we comparethe biasesofthe three estimators discussedintheprevioussection.Asymptotically,theseestimators are all unbiased,but for finite samplesizes thisis not the case. Asstandard MM andtheproposed estimatorbothexist inclosed form, we can analytically study their biases, as we will do via a Taylorapproximation,and findthe conditionsunderwhich the proposedestimatorismoreaccurate.Thereafterweperforma nu-mericalstudytoalsocomparebothestimatorswithML.

4.1. AnalyticalbiascomparisonofstandardMMandtheproposed estimator

ToanalyticallycomparestandardMMandtheproposed estima-tor,weapproximatetheestimatorsbyasecond-orderTaylorseries, asa function ofthe sample mean(

μ

ˆ1) andthe sample variance

(

μ

ˆ2),andthencalculatetheexpectedvalue.Forthederivation,we

referto AppendixA.Here we studytheresults,startingwiththe arrivalrateestimator.ForstandardMM,wefind:

E[ˆ

λ

MM]≈

λ

+

2

λ

n− 1+

2

n.

Fortheproposedestimator,wefind: E[ˆ

λ

NEW]≈

λ

+

exp

(

λ

)

− 1

Firstly,weobservethatbothestimatorshaveapositivebiasforall valuesof

λ

andn.Thisimpliesthatthetruearrivalrateis overes-timatedbybothestimators.Secondly,thebiasesofbothestimators onlydependon

λ

andn,andnoton

μ

c.

Comparing the biases, we find that for

λ

≤ 2.75 the new es-timator is moreaccurate than standardMM, forall samplesizes. For

λ

_{≥ 3.}5,standardMMispreferred.For2_.75_<

λ

_<3_.5,thenew estimator is preferred for smaller samples and standard MM for largersamples.Thesefindingsareinlinewiththeexpectationthat theproposedestimatoroutperformsstandardMMforintermittent demandpatterns.Indeed,since

λ

=2.75correspondstoa6% prob-abilityofzerodemandinaperiod,thenewestimatorispreferred evenforvery‘lightintermittency’.Putdifferently,evenifonly1in 17 periodshaszero demand,thenewestimator still outperforms standard MM. Although both estimators gain accuracy in larger samples, the relative difference remains approximately the same. For

λ

_{= 1}theproposedestimatorisapproximately4timesas accu-rateasstandardMM,whereasfor

λ

=1/4thisfactorhasincreased to17.

We continuewith thedemand size meanestimator. For stan-dardMM,wefind: E[

μ

ˆc,MM]≈

μ

c



1− 1

λ

n



.

Fortheproposedestimator,wefind: E[

μ

ˆc,NEW]≈

μ

c



1−

(

λ

− 2)exp

(

λ

)

+

λ

+2 2

λ

2_n

.

Theseestimatorsbothhaveanegativebiasforallvaluesof

λ

,

μ

c, andn. Thedemand sizemean

μ

c actssolelyasa scaling param-eter. Itfollowsthat for

λ

<2 thenewestimatorismoreaccurate than standardMM.Thescenario with

λ

₌2corresponds toa14% probability ofzero demandina period.We seesimilar dynamics asforthearrival rateestimator.For

λ

=1theproposed estimator isapproximately7timesasaccurateasstandardMM,whereasfor

λ

=1/4thisfactorhasincreasedto170.

There are manyreallife settingswithmoderateto high inter-mittency.Lenguetal.[12] reportthatmorethan80%oftheir stud-ied itemshavean estimated arrivalratebelow 1,themajority of whicharebelow1/4.TurriniandMeissner[13] reportmedian de-mand intervalsof 10 and17 periods, anda maximum of80 pe-riods. In Section 4.2 we present more insights intothe biases of standard MM, ML,and the proposed method by means of a nu-mericalanalysis.

4.2. Numericalstudyofallestimators’biases

Inthissectionweperformanumericalanalysisofthebiasesof standard MM, ML,andthe proposed method,underan exponen-tialcompoundingdistribution.Weselect

λ

=1/16,1/4,and1, rep-resentinghigh,moderate,andlightintermittency,respectively,and presenttheresultsfor

μ

c=2only.Experimentswithother values andtheanalysisinSection4.1 showthat

μ

c actsonlyasascaling parameter. We let the samplesize range fromn=4 to n=200. Each time we draw 1,000,000 period demand histories from the corresponding compound Poissondistribution. Foreach draw, we calculatetheparameterestimatesfor

λ

and

μ

c using(i) the stan-dard MM estimator, (ii)the ML estimator, and (iii)the proposed method.We presentper parametercombination, per samplesize, andperestimator,theaverageofallobtainedestimates.

Two scenarios deserve special attention. The first is the case whereonlyperiodswithoutdemandareobserved.Inpracticeone willnotfitanydemanddistributionthen,butinnumerical exper-imentsitdoesoccasionallyoccur.Allmethodsarethenundefined andtoenableafaircomparison,we setinthiscaseall estimated arrival rate to0 anddo not define demandsize mean estimates.

The other extreme scenario is that with no periods without de-mand.Inthisunlikelyevent,wesettheproposedmethodequalto standardMM.

Fig. 1 showsthe results, with on the left hand side the esti-matedarrivalrates

λ

andontherighthandsidetheestimated de-mand size means

μ

c. Note that in this and all following figures themarkerlocationshavebeenshiftedbetweenestimatorsto en-hancereadability.The scenariowith

λ

=1/16(Fig.1aandb) cor-respondstoa94% probabilityofhavingzerodemandinaperiod. ThestandardMMestimatorconvergesthemostslowlybyfar.For n=200,its arrivalrateestimate isstill 13%toohigh, whereas its demandsize meanestimate is8% toolow. Theproposed method andtheMLestimator performalmost identically,andare consid-erablymoreaccuratethanstandardMM.Forn≥ 50thedifference with the true parameters is below 1%. As the sample size de-creases, the probability of observing no demand at all increases. Especiallyforlowarrivalrates,thisaffectstheresults.

The scenario with

λ

=1/4 (Fig. 1c and d) corresponds to a 78%probability ofhavingzerodemandina period.Standard MM stillshowsasignificant deviationfromthetruevalueforn₌200 (5%for

λ

and2% for

μ

c). The arrival ratesconvergemoreslowly than the demand size means, and the proposed estimator actu-ally(slightly) outperformsML.In Fig.1e andf,thearrival rateis increased to 1, corresponding to a 37% probability that a period haszerodemand.ThestandardMMestimatorstillperformsmuch worse thanthe proposedmethod andML,butthegap issmaller than for lower arrival rates. The proposed method now outper-formsMLbyalargermargin.

Comparing relative accuracies, we find that the ML and pro-posedarrival rateestimator performbest forlow

λ

,whereas the standard MM arrival rate estimator becomes better when

λ

in-creases,although it isstill severely outperformedeven for

λ

=1. Allestimatorsforthedemandsizemeanaremoreaccurateas

λ

in-creases,asfewerperiodswithoutdemandoccur.Directlyrelatedto thisobservationisthequestionoftheoptimalperiodlength(e.g. hour,day,week,ormonth).Differentperiodlengthsleadto differ-entlevels ofintermittency,corresponding to differenttruevalues of

λ

.InSection 6 we studythe resultingeffectoninventory per-formance.

AnexplanationfortheverysimilarperformanceofMLandthe proposedestimatorforlowarrivalrates(hereandforthe Poisson-geometricdistributioninAppendixB)canbefoundintheir defini-tions.TheprobabilitymassofacompoundPoissondistributionfor periodswithnodemandreducestoexp

(

−

λ

)

,whichisminimized at

λ

=0.Forlowarrival ratesmanyperiodshavenodemand.The samplelog-likelihoodfunctionisthesumofthelog-densitiesofall periodsandthereforerelativelyinsensitiveto(few)nonzero obser-vations.Thisinsensitivity alsoholds fortheproposed arrival rate estimator,asthenaturallogarithmisrelatively flataround1(the pointwheren0=n).ThestandardMMarrivalrateestimator,

how-ever, is more sensitive to (even few) nonzero observations, as it usesthesquaredmeanandthevarianceofthesample.Forhigher arrivalrates,thedemandpatternbecomeslessintermittent, bene-fitingstandardMM.ThatisalsowherethebiasofMLdiffersmore fromthatoftheproposedmethod(e.g.Fig.1eandf).

Theresultsofthissectioncanbeusedtovalidatetheaccuracy oftheanalyticalbiasapproximationofSection4.1.Table2 presents thecomparisonforthecase

λ

=1/4,

μ

c=2.Theproposed estima-tor’sexpectation isbetterapproximatedby thesecond-order Tay-lorseriesthan thestandardMM estimator,aswouldbeexpected basedon their differentfunctionalforms. Forn=10,the analyti-cal approximation forstandardMM differs substantially fromthe numericalevaluation.Therearetworeasonsforthis:1)theTaylor seriesapproximationisgenerallymoreaccurateforlargersamples, and2) forsmall samplesizesthe specialcases mentioned above occurmore frequently.Forn≥ 50 thedifferencebetweenthe

an-Fig. 1. Sample Averages of Demand Parameter Estimates, μc = 2 . Vertical Axes Are Truncated for Readability.

alytical andnumericalresultsis small.Theproposed estimatoris well approximated by the Taylor seriesfor all sample sizes. Fur-thermore,underboththeanalyticalandnumericalevaluation,the methodsretainthesameranking.

The proposed method has another important advantage over ML:ithasasimple,closed-formsolution,whereasoptimizationof thelog-likelihoodfunctionforMLinvolvesatwo-dimensional nu-mericalsearch,eachtimeevaluatingafunctionwithinfinitesums. Althoughthisprocedureisnotprohibitivelyexpensive, itdoes re-quiremuchmorecomputationtimethanstandardMMorthe pro-posedestimator.Table3 showstheaveragecomputationtimeper

item(forboththearrivalrateanddemandsizemeanestimate, ex-cludingthesimulationofthedemandsample)forthecase

λ

=1,

μ

c=2(executed on a single Intel Xeon 2.5GHz core).The pro-posedandstandardMMestimatorshowverysimilarandfast com-putationtimesthatmildly increasewiththesamplesize. Contrar-ily,computationoftheMLestimatestakes,forn=10,almost2000 timeslonger than computingthestandard MM orproposed esti-mates,whereasforn₌200ittakesover14000timeslonger.This difference issignificant inlarge-scale applications.Estimating the parametersfor100,000itemseachwith200historicalobservations

Table 2

Analytical and Numerical Expectation Approximations, λ= 0 . 25 , μc = 2 .

Approximated expectation n = 10 n = 50 n = 100 n = 200 ˆ λ, proposed, analytical 0.2642 0.2528 0.2514 0.2507 ˆ λ, proposed, numerical 0.2655 0.2529 0.2514 0.2508 ˆ λ, standard MM, analytical 0.5056 0.3002 0.2751 0.2625 ˆ λ, standard MM, numerical 0.3791 0.2880 0.2709 0.2613 ˆ μc , proposed, analytical 1.9953 1.9991 1.9996 1.9998 ˆ μc , proposed, numerical 1.9701 1.9986 1.9995 1.9997 ˆ μc , standard MM, analytical 1.2000 1.8400 1.9200 1.9600 ˆ μc , standard MM, numerical 1.3784 1.8392 1.9202 1.9596 Table 3

Mean Computation Times per Item for λ= 1 , μc = 2 (Milliseconds).

Computation time (ms)

n = 10 n = 50 n = 100 n = 200 Proposed 0.019 0.021 0.023 0.026 Standard MM 0.019 0.021 0.023 0.025 Maximum Likelihood 3.6 10.7 19.7 36.8

takes61minutesbyML,andlessthan3secondsbytheproposed orstandardMMestimator.

5. Variancecomparison

Inthissectionweusethesamenumericalset-upasinthe pre-vioussection,butreportthevariancesoverthedraws.Fig.2 shows theresults,withthearrivalratesonthelefthandsideandthe de-mand sizemeans on therighthand side.We firstly observethat the proposedmethodandMLperform verysimilarly. Thiscanbe explainedanalogouslytotheirsimilarbiases(seeSection4.2):they differ only intheir treatment ofnon-zero observations, towhich both arerelativelyinsensitivecomparedtothestandard MM esti-mator.SincethevarianceofMLactsasanasymptoticlowerbound anditispracticallyindistinguishablefromtheproposedestimator, we conclude that the new estimator showsnear-optimal perfor-mance.

Forthearrivalrateestimators,standardMMhasasubstantially largervariancethantheothertwomethods.Asthefractionof non-zeroobservationsincreaseswith

λ

,sodoesthedifferencebetween standard MM andthe other two methods. However, for the de-mandsizemeanitcanbeobservedthatonlyfromacertain mini-mumsamplesizeonward,standardMMhasalargervariancethan the other two methods. This reflects a drawback of our numer-ical set-up. If only periods without customer arrivals are drawn, thedemandsizemeanisnotestimated.Thereportedaveragesare therefore biased towards ‘less intermittent’ draws, favoring stan-dard MM.Forlargersamplesizes(where standardMMdoeshave a significantly largervariance thanthe other methods)thisevent virtuallyneveroccurs.

6. Inventorycontrolimplications

For the order-up-to inventory model describedin Section 2.2, we study the service effects when demand parameters are ob-tained according to the three consistent estimators, againfor an exponential compoundingdistribution.Weuseatargetfillrateof 95%andreporttheachievedfillrates.

Fig.3 ashowstheresultsfor

λ

₌1_/16,

μ

c=2,andL=2. Con-sistently withthe results ofSection 4, the ML andthe proposed methodconvergemuchfasterthan standardMM.Theorder-up-to levelthatissetbyanyofthethreemethodsisdrivenbythe over-estimation of

λ

(leading to overshooting the order level and fill rate) andthe underestimation of

μ

c (leading to an undershoot).

For

λ

=1/16,

μ

c=2, andL=2 theunderestimation of

μ

c domi-natesforallthreemethods,leadingtofillratesthatconvergefrom below.Although the overestimationof

λ

ismore severe,this pa-rameteraffects onlytheinventory leveldistribution (viathe lead time demand), whereas

μ

c affects also the distribution of a sin-glecustomer’sdemandsize,whichbothappearinthefillrate ex-pression.WhereastheproposedmethodandMLarealreadywithin 0.1%fromthe95%fillratetarget forasampleof50observations, standardMMonlyattains93.8%evenforn=200.

When

λ

₌1_/4 (see Fig. 3b), the overestimation of

λ

and un-derestimationof

μ

cby theproposedmethodandMLalmost can-celeachother outattheorderlevelcalculation,leadingtoa real-izedfillrateveryclosetoits targetovertheentirerangeof sam-plesizes.InFig.3c(

λ

=1)allmethodsleadtofillrateovershoots, indicatingthat theoverestimationof

λ

isdominant.MLperforms onlymarginally better than theproposed method andthe differ-encesbetweenthethree approachesaresmaller thaninthe pre-viousscenarios.Inthelastscenario(Fig.3d),

λ

isresetto1/4,but theleadtimeisquadrupledto8periods.Comparingthisscenario withthatwhere

λ

₌1_/4andL₌2,weobservethatstandardMM actually performs better under a longer lead time. By extending theleadtime,moreweightisgiventothearrivalrate,sothatthe underestimationof

μ

ciscounteracted.

Summarizingallresults,wefindthatthestandardMM estima-tors show slow convergence of achieved fill rates, especially for highlyintermittentdemandpatterns.MLandtheproposedmethod performclosertotheir targets.Forlow arrivalrates,theproposed method,ML,andstandardMMallconvergetotheoptimalfillrate frombelow,whereasforhigherarrivalratesthispatternshiftstoa convergencefromabove.Theanalysisinthissectionshowsthat– giventhatdemandsarestoredperperiod– thebestperiodlength foroptimalinventory performanceisa non-trivialfunctionofthe arrivalrate,demandsizemean,andleadtime.

7. Robustnessanalysis

In real forecasting and inventory control applications one should notonly estimate demandparameters, butalso selectthe distributional form. In thissection we analyze the robustness of thedifferentmethodsto misspecificationofthedemand size dis-tribution. Wefirst analytically demonstratethe robustnessof the proposed methodand thennumerically compareits performance tothatofstandardMMandML.

Under the demand model specified in Section 2.1, the pro-posedestimatesforthearrivalrateandthedemandsizemeanare ˆ

λ

=− ln

(

n0/n

)

and

μ

ˆc=μˆ_ˆ1

λ.Weexclude,asinSection4.2,the

sce-narioswhere n0=0 orn0=n,astheproposed estimatoris then

undefined. In all other cases,

λ

ˆ and

μ

ˆc are continuous transfor-mationsof n0/nand

μ

ˆ1.Irrespective ofthe (discreteor

continu-ous)compounding distribution, the numberof periods withzero demandoutofnperiodsintotalisbinomiallydistributedwith pa-rametersnandexp

(

−

λ

)

(theprobability thatnocustomerarrives inaperiod).Therefore,E

(

n0/n

)

=exp

(

−

λ

)

.Also,E

(

μ

ˆ1

)

=

λμ

c. Ap-plying the continuous mapping theorem yields that

λ

ˆ→

λ

and

ˆ

μ

c→

μ

c. In conclusion, the proposed estimators have the same closedformandconvergetothetruearrivalrateanddemandsize meanforanycompoundingdistribution.

ThisdoesnotholdforthestandardMMestimator.Forthe pa-rametersofthePoisson-exponentialdistribution,thestandardMM estimatorgivesˆ

λ

=2μˆ21

ˆ

μ2 and

μ

ˆc=

ˆ

μ2

2μˆ₁.Thesecondcentralmoment

of a Poisson-exponential distribution is

μ

2=

λ

(

σ

c2+

μ

2c

)

, where

σ

2

c is the variance of the compounding distribution [40]. There-fore, ˆ

λ

→ 2

(

λμ

c

)

2

λ

(

σ

2 c +

μ

2c

)

.

Fig. 2. Sample Variances of Demand Parameter Estimates, μc = 2 .

Substituting

μ

c=1/λand

σ

c2=1/λ2 fortheexponential distribu-tion, the righthandside simplifies to

λ

,butif thecompounding distributionisnotexponential,thenthisdoesnotholdtrue gener-ally.Similarly,wecanderivethat

ˆ

μ

c→

σ

2

c +

μ

2c 2

μ

c .

Fortheexponentialcompoundingdistribution,therighthandside simplifies to 1/λ≡

μ

c, but again this does not hold true gen-erally. We remark that for other compounding distributions, the

asymptoticproperties ofanyestimator can be derived ina simi-larway.FortheasymptoticanalysisoftheMLestimatorwe refer toÖztürk[30].

The robustness of the proposed estimator extends to com-poundingdistributions withmultipleparameters. Therationale is that it uses P central moments and the probability of zero de-mandforestimatingPparametersofthecompoundingdistribution andthearrivalrate.ThePthcentralmomentofthefullcompound Poissondistributiondependsonly onthefirst P central moments ofthecompoundingdistribution [40].Thisisnotthe caseforthe standard MM estimator, which always uses P+1 moments forP parameters of thecompounding distribution and thearrival rate.

Fig. 3. Achieved Fill Rates (Target Fill Rate 95%). Vertical Axes Are Truncated for Readability.

This makes the MM estimator sensitive to the (nextmoment of the)compoundingdistributionandthereforegenerallyinconsistent underdistributionmisspecification.

Wedemonstratethisforthecasewherethecompounding dis-tribution hastwo parameters, whichwe take to be themean

μ

c andthe variance

σ

2

c.Theproposed estimatoris nowthesolution tothefollowingthreeequations

ˆ

λ

=− ln

(

n0/n

)

ˆ

μ

1=ˆ

λμ

c ˆ

μ

2=ˆ

λ

(

σ

ˆc2+

μ

2c

)

Thefirsttwoequationsdirectlyleadtothesameestimatorsfor

λ

and

μ

c asforaone-parameter compounding distribution. Sub-stitutingtheseintothelastequationgives

ˆ

σ

2 c = ˆ

μ

2 ˆ

λ

− ˆ

μ

2 1 ˆ

λ

2.

Thisclosedformagaindoesnotdependonthespecific compound-ing distribution. Sincethisequation is directlyobtainedfromthe general formula of the second moment of a compound Poisson distribution,andasˆ

λ

→

λ

and

μ

ˆc→

μ

c,thecontinuousmapping theorem impliesthat

σ

ˆ2

c →

σ

c2, irrespective of the compounding distribution. ThestandardMM estimatorinthiscasereplacesthe first equation by the third central moment equation, which de-pendsonthefirstthreecentralmomentsofthecompounding dis-tribution [40]. Therefore, it uses the third moment of the com-poundingdistributiontoestimateits secondmoment.Ifmore

pa-rameters(ormomentsofthecompoundingdistribution) aretobe estimated, then therecursive solution procedureof theproposed estimator(andthereforealsoitsrobustness)readilyextends.

We continue to numericallyanalyze the achieved accuracy of standardMM,ML,andtheproposedestimatorunder misspecifica-tionof thecompounding distribution.In linewiththeremainder ofthepaper,we assumethat thedemandsizesare exponentially distributed. However, we use two different distributions to gen-eratedemandsizes,namelythe(continuous) uniformdistribution withbounds0and2

μ

c,andthe(discrete)logarithmicdistribution (whichwasalsousedbyAxsäter[2])withmean

μ

c.Theproposed estimatorsforthearrivalrateandthedemandsizemeanhavethe sameclosedformforanycompounding distributionandconverge tothetruevaluesby construction.Theotherestimators,however, doshowasymptoticbiasesundermisspecification.Fig.4 showsthe resultsfortheuniformdistribution.

ThestandardMMestimatorfor

λ

showsanasymptoticupward bias of52% inall scenarios, whereas for

μ

c it showsan asymp-toticdownwardbiasof33%inallscenarios.Thismethodtherefore doesnot yield reliable estimates at all when the demand distri-butionis misspecified.MLismorerobust thanstandard MM,but its asymptoticbiasis still considerableforlarger valuesof

λ

.For

λ

= 1/16theerrorisnegligible.For

λ

=1/4werecordforn=200 a0.7%upwardbiasin

λ

anda0.3%downwardbiasin

μ

c.For

λ

=1 thiserrorhasincreasedto+5%for

λ

and-4%for

μ

c.Theproposed method is, in line with our expectation, insensitive to the mis-specifieddemanddistributionandshowsan unabatedaccuracyin comparisontotheanalysisofSection4.2.Alsointhedistribution

Fig. 4. Averages of Demand Parameter Estimates under Distribution Misspecification, Uniform, μc = 2 . Vertical Axes Are Truncated for Readability.

misspecificationcase,thedifferencebetweentheproposedmethod andMLislargerif

λ

isnotverysmall.

Fig.5 showstheresultsforthelogarithmicdistribution.Asthe logarithmic distribution is discrete, it can happen (especially for smallsamplesizes)thatall observationshavethesamevalueand thevarianceiszero.ThestandardMMestimator(forthe exponen-tialcompoundingdistribution)isthenundefined,whichisthe rea-sonwhyinFig.5eandfthebiasofstandardMMisnotdrawnfor n≤ 7.ThisissuedoesnotexistforMLandtheproposedestimator. TheplotsarequalitativelysimilartothoseofFig.4.Asthe logarith-mic distribution’sshape morecloselyresemblesthat ofthe expo-nentialdistribution,theasymptoticbiasesofMLandstandardMM

are smallerin thiscase. Nevertheless,standard MM still showsa largebiasinallscenarios,andfor

λ

₌1wealsoseeaclear asymp-toticbiasinML.Onlytheproposedestimatorisconsistent(andhas thesmallestbiasforallsamplesizes)inallset-ups.

8. Casestudy

In this section we compare the different estimators for com-pound Poisson demand parameters in practice, using a data set of demands for 496 differentcommercial airline parts. This data set containsdemand observations at the individual customer or-derlevel,whichweaggregateto125weeklytotals.Werestrictour

Fig. 5. Averages of Demand Parameter Estimates under Distribution Misspecification, Logarithmic, μc = 2 . Vertical Axes Are Truncated for Readability.

analysis to the 404 parts that havehad at least one demand in the first20weeks.InTable4 wepresentsome descriptive statis-tics.Theestimatedarrivalratesanddemandsizemeansare calcu-latedbasedontheindividuallystoreddemandsofthefull history peritem.Almostallarrivalratesarebelow1,whereasthedemand sizesvarywildly.

Based on theweekly aggregateddemand observations,we set up a comparative analysis forthe arrival rates and demand size means accordingtostandardMM, ML,andtheproposed method. Asabenchmark,wefurthermoreincludethe‘theoreticallybest’ es-timatesbasedontheindividualorders,which,ofcourse,cannotbe exploitedusingaggregateddata.Lenguetal.[12] establishedthata

compoundPoissonprocess withgeometric compounding distribu-tionhasagoodfittothisdata.Weadoptthisdistribution(whichis furtherstudiedinAppendixB),butinlinewiththerestofthis pa-per,we alsotest theexponentialcompoundingdistribution. Start-ing withthe first 20 weeks and every time increasing the time horizon by one week, we estimate for every horizon the arrival rateanddemandsizemeanforallitems.WereporttheMean Ab-solute Percentage Error (MAPE) betweentheseestimates andthe ‘theoreticallybest’estimatorovertheentiresample(foreverytime horizon,asanaverageacrossallitems).Theproposedmethodand the‘theoreticallybest’methoddonotdependonthespecific com-pounding distribution. Furthermore, the ML results do not vary

Table 4

Case Study Data Set Descriptive Statistics. No. Items No. Weeks Estimated Arrival Rate

404 125 Min Max Mean Median ≤ 0 . 25 ≤ 0 . 5 ≤ 1 0.04 2.87 0.29 0.19 65% 86% 97% Estimated (Individual) Demand Size Mean

Min Max Mean Median ≤ 2 ≤ 10 ≤ 100 1.07 802.5 21.42 5.34 6% 74% 94%

Fig. 6. Mean Absolute Percentage Errors Between the Different Methods and the Theoretically Best Estimator for Case Study Data.

visibly betweenthe exponentialandgeometric compounding dis-tribution. Therefore, only for standard MM do we report the re-sultsunderboth distributions separately.Theresultsare depicted inFig.6.

By construction, the MAPEofthe individual demands estima-torconvergesto 0andtheevolutionofthisestimatorthroughout the samplebenchmarkstheother methods.ThestandardMM es-timator clearlyperformsworse whentheexponentialdistribution isassumed.However,inbothcasesitperformsconsiderablyworse thantheothermethods.Forthearrivalrate,theproposed estima-torslightlyoutperformsMLby0.2%forthefullsample.Forsmaller samples, the differences betweenthe estimators are smaller, but their performance rankingdoesnot change. Forthe demandsize mean,standardMMevendeteriorateswhenthesamplesizepasses 90 weeks.The theoretical resultsfromtheprecedingsectionsare confirmedalsointhispracticalcasestudy:standardMMperforms byfarworstandismostsensitivetotheexactspecificationofthe compounding distribution,andthe proposedestimatormarginally outperforms ML. It is furthermore established that the accuracy loss due to storing only weekly data increases with the sample size.Ifall 125weeksareused,then theproposedestimatorloses approximately10%inaccuracycomparedtousingindividual trans-action data. Standard MM,contrarily, loses approximately 35% or 40%,dependingonthecompoundingdistribution.

9. Conclusions

The inventory control literature provides hardly anyguidance on obtaining demand parameters from period demand data, al-though companies typically do store historical demand observa-tions per period andsoftware packages also use period demand as input for their forecasts. Period demand forecasts cannot be directly transformedintocompoundPoissonparameter estimates, eveniftheforecastingprocedureyieldsseparateestimatesforthe time betweenperiodswithpositivedemands andtheaverage

pe-riod demand size. Misusing these to obtain demand parameters leads – also asymptotically – to severely biased estimates, over-shotfillrates,andthusexcessiveinventoriesandassociatedcosts. Itisimportanttoobservethatthisflawoccursifanyperiod fore-castingprocedureisusedtoestimateindividualcustomerdemand parameters. We have addressed this lack of connection between the forecastingandinventory control literature, andperformed a comparativestudyofdifferent compoundPoissonparameter esti-mators.WeproposedanalternativeMMestimatorthatusesa non-standardchoice ofmoments, reflectingtheintermittent natureof thedemand.Thefractionofperiodswithoutdemandisusedto es-timatethearrivalrate,andtheestimatedarrival rateiscombined withfurthercentralmoments ofthe(full)compoundPoisson dis-tributiontoestimatetheparametersofthecompounding distribu-tion(thedistributionoftheindividualdemandsizes).

OurresultsconfirmedthatthestandardMMestimator(the de-fault option in inventory control) has severe biases for a wide rangeofintermittentdemand patternsandsamplesizes,andhas a larger variance than the other studied methods. The proposed method’s deviation from the true value of both parameters is smallerovertheentirerangeofintermittentdemandpatternsand samplesizesconsidered. In severalscenarios (particularly also in ourcasestudy)theproposedestimatorhasalowerbiasthanML, whereas in the other scenarios it performs almost indistinguish-ably. A key advantage of the proposed methodis 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 neitherstandard MM, nor MLpossesses)isparticularlyrelevantinpractice.Theclosed-form natureoftheestimatorisusefulinlarge-scaleapplicationswhere manyitemsarecontrolledsimultaneously andcomputationspeed is an issue.Both practitioners and applied researchersshould be aware of 1) the accuracy loss dueto storing demand only peri-odically,and2)theconsequencesoftheirestimatorchoiceon ac-curacyandinventoryperformance. Itis furthermoreimportantto