Stock control analytics: a data-driven approach to compute the fill rate considering undershoots

One of the most frequently used inventory policies is the order-point, order-up-to-level (s, S) system. In this system, the inventory is continuously reviewed and a replenishment request is placed whenever the inventory position drops to or below the order point, s. The variable replenishment order quantity and the variable replenishment cycle characterize the system by the use of complex mathematical computations. Different methodological approaches diminish the mathematical complexity by neglecting the undershoots, i.e., the quantity that the inventory position is below the order point when it is reached. In this paper, we conceptually and empirically analyse the bias that neglecting the undershoots introduces into the estimation of the fill rate. After that, we suggest a new methodology developed under a data-driven perspective that uses a state-dependent parameter algorithm to correct such a bias. As a result, we propose two new methods, one parametric and the other nonparametric, to enhance the fill rate estimate. Both methods, named analytics fill rate methods, remove the bias that neglecting the undershoots introduces and are used to illustrate the practical implications of this hypothesis on the performance and design of the (s, S) system. This research is developed in a lost sales context with simulated stochastic and i.i.d. discrete demands as well as actual sales data.


Introduction
Inventory management plays a critical role in the success of most enterprises because of its direct implication for decision-making both at the operational and strategic levels. In the current markets, with more demanding customers and higher costs, effective inventory management systems are essential. Therefore, developing an intelligent inventory management system and transforming inventory data into valuable insights becomes a key factor in achieving competitive advantages (Stefanovic 2015;Zhou et al. 2017).
The core purpose of inventory management is to address two important key issues: (1) when to launch a replenishment order, and (2) the size of this replenishment order. The control parameters of inventory policies are devoted to determining these issues according to a specific objective, normally using a cost or a customer service criterion. Even though the cost approach often leads to less mathematical complexity, the use of the service criterion avoids the difficulty of accurately estimating the stockout cost due to the presence of intangible components such as the loss of customer goodwill and the diminishing of future business (Gutgutia and Jha 2018) opportunities. For that reason, this paper follows a customer service approach that introduces a control parameter known as the service level that becomes a constraint in determining the control parameters of the system (Silver et al. 2017;Escalona et al. 2019). The service level addressed in this paper is the fraction of the demand that is immediately satisfied from the stock without causing a shortage (Brown 1962), which is known as the fill rate. Despite its simple definition, its calculation hides complex mathematical nuances, especially when the system is managed by a continuous review system (the status of the inventory is known at any given moment) and unfilled demand is lost, as happens in highly competitive sectors such as retail, machinery spare parts, service sector, or e-commerce of fixed-date services (Gruen et al. 2002;Breugelmans et al. 2006;Diels and Wiebach 2011;Disney et al. 2021).
This paper focuses on the order-point, order-up-to-level (s, S) system due to its versatility and extended use in practice (Caplin and Leahy 2010;Bijvank et al. 2015). By definition, this system assumes continuos review and, in a stochastic demand context, both the replenishment cycle, i.e., time elapsed between two consecutive order deliveries, and the order quantity, are variable; therefore, simultaneous determination of the parameters is necessary to guarantee the mathematical optimality of the policy once a service objective is set. However, as Silver et al. (2017) point out, in practice, the most common approach consists of assuming that all the demand transactions are unit sized (see for example Silver 1970;Vincent 1983;Platt et al. 1997;Agrawal and Seshadri 2000;Axsater 2000;Larsen and Thorstenson 2014). This means that the inventory position always exactly reaches the order point, s (known also as reorder point or ROP), and therefore, the order quantity is constant and equal to the order-up-to-level, S, minus the order point, s, (i.e. S − s). However, in real situations, the inventory position may not exactly be at the ROP but at a certain amount below it, known as the undershoot at ROP. Neglecting undershoots greatly reduces the mathematical complexity of the problem but can severely reduce the service level of the system (Schneider 1981), which leads to increased stockout costs, especially when managing lumpy or erratic demand items for which the unitary demand assumption is not appropriate (Cardós and Babiloni 2011a).
In the research literature, we find several studies that address the complexity of estimating undershoots. We find quite a number of works that calculate the undershoots based on the asymptotic results of renewal theory (Schneider 1981;Tijms and Groenevelt 1984;Baganha et al. 1999;Johansen and Hill 2000;Kouki et al. 2009). Other authors, in contrast, compute probability distribution functions to estimate the expected value of the undershoots and incorporate it into service measures (Cohen et al. 1988;Baganha et al. 1996;Gutierrez and Rivera 2021). However, in general, these studies approximate the value of the undershoots, which may have an impact on many aspects of inventory management. Regarding the fill rate estimations in service models, Schneider (1981), Tijms and Groenevelt (1984), Moors and Strijbosch (2002), Silver et al. (2009) andSilver et al. (2012) suggest different approximations of the fill rate that consider undershoots under the backordering assumption. In the lost sales context, only Cohen et al. (1988) present a fill rate method in a order-point, order-up-to-level (s, S) system; however in this approach, the undershoot estimations are calculated from an approximate distribution function that may bias the fill rate estimation.
Then, although the (s, S) system is widely used in practice, there is no optimal approach to establish the control parameters when the fill rate is the service constraint of the system and the unfilled demand is lost. The objective of this research is to propose a new methodology to obtain the fill rate based on a bias analysis and a correction of the undershoot assumption. The methodology follows a data-driven perspective that avoids any assumption regarding the undershoot distribution and is suitable for service models.
The implementation of business intelligence affects the agile efficiency of the supply chain (Kaur 2021). In general, the digitalization of companies offers a large amount and variety of data (Big Data) that transforms the way supply chain management works (Wang et al. 2016). In particular, Gandomi and Haider (2015) indicate the need to develop appropriate and efficient methods to leverage such data, and they review big data analytics areas such as text analytics or predictive analytics. In this sense, this work is an example of how to extract more value from inventory data and proposes use of a novel area called stock control analytics, where data-driven methods, which include machine learning and artificial intelligence techniques, can enhance the state-of-the art inventory management methods.
Although the nomenclature may be new, the use of data analytics tools to overcome the limitations associated with theoretical inventory models is not that new. For instance, Strijbosch et al. (2000) correct the bias introduced by using demand forecasts instead of statistical demand distribution first moments, which are unknown prior to the commencement of any inventory processes. That article proposes a heuristic method based on simulations, which substantially reduces the bias. Novel methods based on this data-driven perspective have been recently published to deal with the uncertainty of demand distributions. Huber et al. (2019) revisit the newsvendor problem subject to a target cycle service level by employing machine learning and quantile regression approaches that circumvent the need to assume a statistical distribution for the demand. Regarding the uncertainty of statistical demand distribution, Trapero et al. (2019a) investigate parametric and nonparametric methods to enhance the safety stock estimation problem for a certain target cycle service level. Their results show that nonparametric approaches such as kernel density estimators provide a good performance in terms of cost and service for lower lead time values, and parametric approaches such as GARCH outperform the rest of the methods for higher lead time values. Trapero et al. (2019b) continue with this line of research by proposing a combination scheme of GARCH and Kernel whose weight parameters are determined by minimizing the tick-loss (newsvendor) linear asymmetric cost function. Additionally, regarding the performance of machine learning forecasting methods versus traditional counterparts, Spiliotis et al. (2020) compared both methods for daily SKU demand forecasting. The work highlights the improvements possible with machine learning methods.
Previous articles show how data-driven approaches offer versatile solutions to address the limitations of traditional theoretical assumptions. In this work, following the same philosophy, we focus on the assumption of neglecting the undershoots on the (s, S) system when the demand is modelled by any discrete distributions for a lost sales case. To the best of the authors' knowledge, this particular topic has been overlooked in the literature. Basically, this research aims to determine the extent of the bias introduced by such an assumption and proposes a methodology to reduce it. The data-driven approach employed to cope with this problem is the state-dependent parameter (SDP) estimation technique. SDP estimation involves the nonparametric identification of the state dependency using recursive methods for time variable parameter estimation, which allow for rapid (state dependent) parametric change (Young et al. 2001). SDP estimation has been successfully applied in a supply chain context in Trapero et al. (2011), particularly to correct judgemental forecast biases. Once the nonparametric estimation between the state and the parameter is available, another advantage of the SDP procedure is that it reveals any potential nonlinear pattern between the state and the parameter and can thus be further parametrized, obtaining a fully self-contained model with just a few parameters (Young et al. 2001). Following this methodology, we propose two methods to estimate the fill rate, one parametric and another nonparametric; the methods remove the bias and outperform the initial approach. These methods are named analytics fill rate methods, by the authors, to distinguish them from the conventional formulas.
The rest of this paper is organized as follows. Section 2 describes the problem the paper addresses and states the main assumptions and the notations. Section 3 introduces the initial approach of the fill rate formula, which neglects the presence of undershoots. This section provides insight of the potential bias introduced by neglecting the undershoots. In Sect. 3, the experimental setup is also detailed. Section 4 explores the SDP approach to correct the bias associated with the initial fill rate estimation and proposes the new analytics fill rate methods. Section 5 assesses the methodology using a case study. Practical implications and a discussion of the results are presented in Sect. 6. Finally, the main conclusions are drawn in Sect. 7.

System description, notation, and assumptions
This paper considers a single-echelon, single-item inventory system where the demand is stochastic. The system is modelled by a discrete time model where the demand arriving during a given interval of time comes from a discrete distribution. The stock is controlled following a continuous review order-point, order-up-to-level (s, S) system for the lost sales case. In this system, when the inventory position is at or below the ROP, a sufficient amount is ordered to raise the inventory position up to the order-up-to-level, S. The replenishment order is received L units of time after being launched. We consider that the replenishment cycle (also referred to as just, the cycle, further on) is the time between two consecutive deliveries and starts at order delivery. Figure 1 shows an example of the evolution of the on-hand stock and the inventory position when the system is not out of stock (a) and when it is out of stock (b).
The notations in Fig Q t = order quantity at instant t (units), β A = achieved fill rate, β 0 = initial fill rate, β i = estimated fill rate for method i, i = LR, SDP.
The general assumptions of this paper are as follows: (1) the lead time, L, is constant and known; (2) there is never more than one outstanding order, which implies that the condition s < S − s is always true; (3) the demand process is considered stationary and i.i.d., and is defined by any discrete distribution function; (4) unfilled demand is lost and; (5) z t is the on-hand used to satisfy d t . Note that assumption (3) is generally accepted in inventory research since, if it is not the case, the numerical difficulties are insurmountable (Schneider 1981), especially for the lost sales case (see, for example, Hadley and Whitin 1963;Cohen et al. 1988;Johansen 2001).
This research considers that time is discrete and is organized into a numerable and infinite succession of equally spaced instants, which is tantamount to assuming that each transaction is updated in a specific instant of time considering a sequence of discrete and infinitely small intervals. This procedure is indeed the way IT systems that are used to control continuous inventories work, i.e., they report every transaction that occurs in a specific instant of time.

Problem description
The (s, S) system implicitly states that both the replenishment cycle and the order quantity are variable. These two characteristics make it truly difficult to find a mathematical approach with the potential for a practical implementation.
First, the size of the replenishment order depends on the stock level when the ROP is reached, i.e., Q = S − z. In a stochastic context, to compute the expected value of z τ , bearing in mind that the evolution of the on-hand stock can be modelled as an ergodic Markov chain, it is necessary to know the probability transition matrix of the on-hand stock between two consecutive replenishment cycles M in such a way that lim n→∞ P(z) ⋅ M n = P(z 0) , where P(z) is an arbitrary vector and P z 0 is the principal left eigenvector of M , i.e., the probability vector of the on-hand stock levels at order delivery. This procedure is very time-consuming, especially for large values of S, since it determines the dimension of the transition matrices (Cardós et al. 2017). If the undershoots are negligible such that z τ = s, the order quantity is Second, the estimation of the total demand and concretely D is not straightforward and depends on the expected values of the stock levels and on the unknown number of periods until the ROP is reached, τ. To delve into this problem, please refer to Cardós and Babiloni (2011b). As a result, the replenishment cycle is variable, which, in a stochastic demand context, hinders the calculation of the demand distribution from the order delivery until the ROP is reached. Nevertheless, if the undershoots are neglected, the demand consumed from the beginning of the cycle until the order point is reached is the difference between the on-hand stock at the beginning of the cycle and the ROP.
For these reasons, in practice, "the values of the control parameters are set in a rather arbitrary fashion" (Silver et al. 2017), and the common derivation of the (s, S) system consists of assuming that all the demand transactions are unit sized or, stated another way, that undershoots at ROP are not possible.
However, what does neglecting undershoots mean? The answer is related to the expected value of the stockouts (also known as expected shortages) during the replenishment cycle. The assumption that the replenishment order is launched at the precise moment the ROP is reached implies that the on-hand stock at launch may be greater than it actually is, since z τ ≤ s, in such a way that there is a higher probability of stockouts than expected. In the following sections, we describe the bias that neglecting undershoots introduces into the estimation of the fill rate as a consequence of underestimating the stockout probability of the system; we then propose a data-driven approach that eliminates this bias. The results are validated using a simulation since there is not an exact method able to compute the fill rate in the context of our research.

Conceptual analysis of the bias
The fill rate is commonly defined as the fraction of demand that is immediately fulfilled from the on-hand stock. However, this simple-looking definition hides many technical details and nuances that are often overlooked in the literature. In fact, the fill rate has been simplified through the traditional approximation which computes the fill rate in terms of units short, i.e., as the complement of the ratio between the expected shortage per replenishment cycle (ESPRC) and the total expected demand per replenishment cycle (EDPRC) Vis 2011, 2012;Guijarro et al. 2012; Babiloni and Guijarro 2020), i.e., Although the definition of the fill rate does not introduce any assumption regarding the inventory system, as pointed out before, one possible approach to avoid the computational complexity that undershoots introduce is to assume that they are negligible at ROP. This assumption implies that, (1) the order quantity is constant and, (2) the system will only be out of stock during the lead time. Figure 2 shows the evolution of the stock levels when the undershoots are neglected, where both features can be easily appreciated.
When the stock exactly reaches the ROP, the system is only out of stock if D L > s, and therefore, the stockouts only take place during the lead time. Hence, the expected shortage per replenishment cycle is straightforwardly computed as To compute the expected total demand per replenishment cycle (EDPRC), we sum the accumulated demand from the beginning of the cycle until the ROP is reached, D, and the accumulated demand over the lead time, D L . The latter is simply computed when the demand distribution is known. However, to compute D τ , we need to know the on-hand stock balance at the beginning of the cycle, i.e., at order delivery, which, in a lost sales context, is obtained by z 0 = S − s + E[s − D L ] + . Therefore, taking into account that D τ must guarantee that the ROP is exactly reached, then Consequently, the expression to estimate the fill rate when the undershoots are neglected (named initial fill rate, β 0 , further on) for the lost sales case that applies to any discrete distribution is: It is clear that conceptually neglecting undershoots leads to a systematic overestimation of the fill rate and, therefore, introduces a bias in its estimation, which can be clearly appreciated in the following example. If the stock at t is one unit above the ROP, and demand at t + 1 is 2 units, the ROP is reached, and a replenishment order is placed. The real stock at launch is s − 1 units. Therefore, the on-hand stock remaining on the shelves available to meet the demand during the lead time is one unit less than what is assumed if the undershoot is neglected. Thus, the initial fill rate is higher than the real rate, and thus, the stockouts are larger than expected. This simple example shows why neglecting undershoots introduces a significant bias in the estimation of the fill rate. (1)

Empirical analysis of the bias
To analyse the bias of the initial fill rate, we designed an experiment based on simulations as follows (see Fig. 3): (1) we combine a set of input parameters that define the inventory policy and the simulated demand; (2) for each combination of the input data, we simulate a (s, S) inventory system using the Monte Carlo method and calculate the achieved fill rate; (3) for each input data combination, we also compute the initial fill rate using expression (1). Note that inventory simulation is an adequate approach given the complexities of the system (Chinello et al. 2020).
With regard to the input data, an extensive range of values for s, S and L is selected to provide realistic values of the fill rate (from 0.5 to 0.99). Regarding the demand, we select a set of parameters (r, θ) for the negative binomial distribution based on the consideration of four demand categories suggested by Syntetos et al. (2005): smooth, lumpy, intermittent and erratic, where r is the number of successes and θ is the probability of success. Moreover, the negative binomial distribution reproduces some characteristics typically seen in demand data; for instance, it is positive and asymmetric, which models potential marketing campaigns. Tables 1 and 2 present the set of data used in the experiment, which considers every feasible combination of these values per factor, that leads to 1360 different cases (excluding cases where s ≥ S − s to ensure that, as we assume, there is never more than one  Order-up-to-level S = 5, 7, 12, 15 Lead time L = 1, 2, 3, 4  outstanding order). The inventory system simulation follows the diagram (Fig. 4), and for each input data combination, the fill rate achieved (β A ) is obtained as the average fraction of the complement of the unfulfilled demand over the total demand in every replenishment cycle when considering 20,000 consecutive periods, i.e.: where N indicates the total number of replenishment cycles. To assure the consistency of the results, Monte Carlo simulations undertake thirty replications of each case using the average of these replications as the final achieved fill rate. Figure 5 compares the achieved (β A ) and the initial (β 0 ) fill rates plotted as red circles and illustrates the bias introduced in the calculation of the initial fill rate due to neglecting the undershoots. The larger the distance between the red circles and the blue line, the larger the bias. Note that the initial fill rate always overestimates the achieved fill rate, as explained in Sect. 3.1. What is empirically found is that the bias size is not constant, and it seems to increase over the central part of the figure.
Therefore, the problem is to find the bias between β 0 and β A , which appears to have a nonlinear pattern.

Analytics fill rate methods
To model the bias shown in the previous section, we assume that the data are crosssectional, i.e., the relationship between the achieved fill rate and the initial fill rate does not depend on time. The initial number of experiments was 1360. We focused on fill rates greater than 50%, filtering the number of experiments to 1011. The training and test sets are divided by random sampling, where 70% of the data (708 experiments) are used for estimation purposes and the rest of the data (303 experiments) are used for validation. Figure 5 shows that the bias between the fill rates is nonlinear and is greater for fill rates between approximately 0.55 and 0.9, where the maximum bias is reached at approximately 0.6 and 0.75. To model that bias, this work adopts an SDP approach. (

Fig. 4 Diagram of the simulation
The approach requires that we define the state and the parameter that results in a variance of the state. An initial model could be as follows: where α 1 is the state that somehow depends on the parameter β 0 . Note that we do not exactly know the relationship between the state (α 1 ) and the parameter (β 0 ), and it will be the SDP algorithm, using the fill rate simulated data, that will provide a curve to give us some insights about the pattern of such a dependency. In this sense, SDP approximates β A .

A nonparametric approach: the SDP algorithm
The SDP approach can be seen as an extension of the time variable parameter estimation (TVP) (Young 2012), where the unknown parameters are slowly variable with time. Typically, TVP estimations are represented by a two-dimensional state (level and slope) to model its stochastic behaviour by means of a generalized random walk (GRW), which is a generic model to unify the different versions of the random walk as an integrated random walk (IRW), random walk (RW) and smoothed random walk (SRW). More information about the use of TVPs in time series analysis and forecasting can be found in Young et al. (2001). The SDP approach employs a similar procedure but the main difference is that the parameters evolve stochastically with respect to another variable instead of time. Such parameters are denoted as state-dependent parameters (SDPs). This can be achieved by "sorting" the data in a nontemporal order. If the new ordering provides a smoother and less rapid variation in the SDP, a GRW will be able (3) SDP = 1 0 ⋅ 0

Fig. 5 Achieved fill rate versus initial fill rate (red circles)
to describe its evolution in the transformed observation space. Utilizing a fixed interval smoothing (FIS) estimation, the estimated SDP can be "unsorted" to the original order.
In our particular case, the SDP α 1 can be sorted with respect to β 0 . This new organization of data is indexed by k, and thus, the stochastic dynamics of the SDP can be expressed as an integrated random walk defined in a state space framework as follows: where α 1 is the level and * 1 is the slope. Small variations with respect to the new sorting are introduced with the random Gaussian noise * (k) with zero mean and the constant variance σ 2 α . Such a variance is often referred to as a hyperparameter to differentiate it from the states that are the main object of the estimation analysis. The observation equation needed to complete the state space representation is: The last term is the typical error component that assumes a Gaussian distribution with zero mean and constant variance σ 2 .
Fortunately, these routines are already implemented in the CAPTAIN toolbox of MATLAB available in http:// capta intoo lbox. co. uk. See "Appendix" for more information about the recursive filtering/smoothing algorithms used. Figure 6 shows the estimated SDP α 1 (β 0 ) sorted with respect to the parameter β 0 for the training set. This figure shows the nonlinear pattern of the SDP coherently with the nonlinear bias found between β A and β 0 in Fig. 6. (4) 1 (k + 1) * 1 (k + 1)  Figure 7 depicts the SDP estimation of the fill rate (β SDP ) together with the initial and achieved fill rates for the training set. That figure clearly shows how the SDP approach is capable of considerably reducing the bias. Essentially, the fill rate estimates provided by β SDP estimate the mean of β A , correcting the systematic deviation found for its initial counterpart β 0 .

Parameterization based on the SDP graph
The SDP approach provides a nonparametric estimation of the relationship of α 1 (β 0 ) based on a graph, as shown in Fig. 6. Furthermore, this approach can be complemented by parameterizing such a curve and proposing another parametric alternative. Figure 8 shows two polynomial (quadratic and cubic) functions that could be used to fit the resulting nonlinear SDP α 1 (β 0 ) curve.
Considering that, the following model can be proposed: where i for i = 0, 1, 2, 3,4 are constants that can be estimated by means of a linear regression and is the typical error component that assumes a Gaussian distribution with a zero mean and a constant variance. Again, LR approximates A . It is important to note that, despite the nonlinearities involved in the initial fill rate bias, thanks to the SDP identification, we have arrived at a nonlinear model that is linear with respect to the parameters i . Note that such parameterization makes the model fully self-contained and clearly reveals the nature of the graphically identified nonlinearity with just a few parameters.
Fill rate estimations provided by the initial ( 0 ) and the SDP ( SDP ) approaches for the training set Table 3 summarizes the estimation results for the quadratic and cubic polynomials. The table shows that all parameters are statistically significant for a significance level of 1%. Although both approaches perform well, the cubic polynomial achieves a smaller RMSE in the estimation and a slightly better goodness of fit ( R 2 ).
To validate that the parameters can be used to correct the bias with data that have not been used in the training set, we used those estimated parameters in the data test set. Figure 9 shows the fill rate estimates provided by the nonparametric SDP (β SDP ) and the parametric cubic polynomial model (β LR ) compared to the achieved and initial fill rates. That figure shows that the performance provided by β SDP and β LR is very similar, corroborating the utility of both models.

Error assessment
We can define the fill rate error as: where β ij for i = 1, 2, 3 is the fill rate estimated by the initial equation, the nonparametric SDP and the parametric linear regression, respectively. The index j = 1, 2, …, J is the number of fill rate observations for the entire test set J. The performance of the different fill rate estimates can be summarized by the mean error (ME) and the root mean squared error (RMSE), respectively, such that: Note that the ME measures the error bias, whereas the RMSE focuses on the error size. Figure 10 shows the ME and RMSE obtained for the different fill rate estimation techniques. First, the high positive value of ME for β 0 quantifies the overestimation produced by the expression (1). Note that such a formula also displays a high error variability as measured by the RMSE. This figure also clearly shows how the parametric approach (β LR ) and the SDP (β SDP ) alternatives outperform the initial (β 0 ) in Fig. 9 Fill rate estimations provided by the initial formula (β 0 ), SDP approach (β SDP ) and parametric model (β LR ) for the test set terms of the bias (ME) and error variability size (RMSE). On the other hand, the differences found between β SDP and β LR are minimal.

Case study
To verify the performance of the proposed methods for a real dataset, sales data previously used by Barrow and Kourentzes (2016) were utilized as the demand estimates. This dataset belongs to a major UK fast-moving consumer goods manufacturer specializing in the production of household and personal care products. In total, 229 products at the SKU level were available, with 173 weekly observations per product. According to Barrow and Kourentzes (2016), there are no seasonal SKUs, and only a minority (21%) exhibit a slightl trend. To calculate the achieved fill rate for real data, we introduced the demand data of every SKU in our simulation, indicating the values of S, s and L for each product, and following the procedure summarized in Fig. 4, we calculated the fill rate with expression (2).
Proceeding as in the previous simulations, 70% of the data are held as a training set, and 30% are held as a test set. Figure 11 shows the initial (β 0 ) and achieved (β A ) fill rates for the training set of the case study. As previously occurred with the simulated demands, there is a nonlinear positive bias of β 0 with respect to β A . That bias is higher around fill rates of 0.8 and 0.9. Figure 12 depicts the results obtained in the test set, where β SDP and β LR have been reestimated using the training dataset. Note that β LR was computed by means of a quadratic polynomial function. In general, the bias of β SDP and β LR was significantly  Relationship between the initial fill rate (β 0 ) and achieved fill rate (β A ) for the training set of the case study data Fig. 12 Fill rate estimations provided by the initial formula (β 0 ), SDP approach (β SDP ) and parametric model (β LR ) for the test set of the case study dataset reduced. The bias reduction is quantified in Fig. 13. Again, the ME and RMSE of β 0 with regard to the proposed alternatives are larger.

Discussion and practical implications
The previous sections showed that analytics fill rate methods successfully remove the bias associated with the initial fill rate formula. In this section, we will illustrate the managerial implications and benefits of using them in practice. Essentially, we will describe two possible applications: (1) a system performance measurement and (2) a policy design. In both examples, we will use the parametric alternative (β LR ), since it provides similar results when compared with the nonparametric counterpart (β SDP ).

System performance measurement
For companies, one possible application is to measure the performance of a determined (s, S) pair of parameters in terms of the fill rate. Consider a manager that sets a fixed value of S for physical restrictions, and s is defined judgementally. The manager would like to know beforehand what the resulting fill rate for those chosen parameters will be, given past sales as demand estimates and a given lead time. By using the initial formula, they would obtain 0 . However, as previously demonstrated, such a formula overestimates the real fill rate, and the system will be less protected than the manager expects with an unexpected increase in stockout costs. Thus, it should be corrected by means of LR .
An example of that application is shown in Fig. 14. Given the sales data from SKU 126, S = 1000 and L = 2, the figure shows that for different values of s, the 0 overestimates the fill rate achieved A . In contrast, LR follows A in a closer fashion.

Parameter design
Another application of the proposed methodology is to design an inventory policy with the fill rate as a constraint. In other words, given a target fill rate, which values of s and S should be chosen? In this case, we have two options. First, if we do not apply the bias correction, then 0 = target , and by means of Eq. (1), s can be obtained, denoted by s 0 . Alternatively, to avoid the bias associated with 0 , we set LR = target , and we can obtain ′ 0 from expression (9). Note that, unlike Eq. (6), in the expression (9), LR acts as an input, and ′ 0 is the output. Proceeding in the manner shown in Sect. 4, ′ i for i = 0, 1, 2, 3,4 should be reestimated with the training dataset. Note that the parameters in (6) and (9) are not the same.
With the new value of ′ 0 together with Eq.
(1), we can obtain the value of s, which is denoted by s LR . Next, these values of s 0 and s LR are fed into the simulation to obtain A (s 0 , S) and A (s LR , S) , whose values are subject to a comparison with respect to the initially defined target .  Table 4 shows various examples for three different SKUs of the real dataset. For instance, the first row of the table is an example for SKU 126, given L = 2, S = 1000 and target = 0.75. The calculation of the reorder points is s 0 = 104 and s LR = 300, depending on whether the bias is corrected. These ROPs are fed as inputs into the simulation, and the fill rates achieved are 0.72 and 0.76, respectively. Since the fill rate achieved by A s LR , S is closer to target , it corroborates the usefulness of the proposed method to more effectively design the inventory policy, which means a reduction in stockout costs. Similar conclusions can be drawn regarding the rest of the SKUs in the table.

Summary and conclusions
This paper focuses on the versatile order-point, order-up-to-level (s, S) system in the lost sales context. The control procedure of this system consists of continuously reviewing the status of the inventory in order to know exactly when the inventory position is at or below the ROP in order to launch a replenishment order to raise the inventory position to the order-up-to-level, S. Despite this system being widely used in practice, it is normally implemented assuming that the undershoots at ROP are negligible, i.e., assuming the inventory position always exactly reaches the ROP, which is only possible if demand transactions are unit sized. However, if demand is not unitary, the assumption regarding neglecting the undershoots leads to a bias in the performance of the system. This paper shows both conceptually and experimentally that the initial formula to compute the fill rate is biased. It can provide an overestimate of almost 7% in the estimation of the fill rate and lead to wrong decisions that directly affect the performance and costs of the inventory system. This paper proposes a new methodology from a data-driven perspective that uses a nonparametric SDP algorithm to model the bias and subsequently suggests two analytics fill rate methods that outperform the initial approach. Note that, despite the unknown complex nonlinearities involved in the initial fill rate bias, the SDP has revealed the nature of such a nonlinearity, which has been subsequently parameterized by a parametric nonlinear model, linear in the parameters, that allows the use of linear regression. The importance of the analytics fill rate methods is that both are unbiased and easily implemented in practice. Although the SDP approach can be seen as a data-driven tool, it is important to distinguish it from other black-box data-driven tools typically associated with machine learning alternatives. In this work, the SDP approach was employed to learn the nature of the nonlinearities involved in the initial fill rate bias by neglecting the undershoots in the (s, S) system. Once those nonlinearities were identified by means of a curve, the expert model could develop a parametric solution taking advantage of the solid theory that underpins parametric theory. In this case, a simple linear regression was enough to model a complex nonlinearity. Note that this philosophy follows the Data Based Mechanistic approach proposed by Young (2012), where the author explains that "a model should not just explain the time series data well, but it should also provide a mechanistic description of the system under investigation".
Practitioners interested in applying the methodology proposed should proceed as follows. First, they should compute the potential bias between their theoretical fill rate expression used and their achieved fill rate; if that bias is similar to the one presented here in the simulations and case study, then they can directly utilize the parametric model ( LR ) . If the bias found does not follow a similar nonlinear pattern, then the previous SDP identification step should be carried out. Note that the improvements obtained per SKU might mean a significant reduction in the stockout costs when applied to a large portfolio of SKUs.
Further research should analyse the effect on replacing the statistical demand distributions in the fill rate initial formula with the probabilistic forecasts of real demands. In addition, further research is also needed to extend these promising results to other systems.

Appendix
This appendix shows the recursive filtering/smoothing algorithms that are usually employed to reflect the TVP evolution (Young 2011) and how they have been extended and modified for SDP estimation. More details about this methodology can be found in Young (2011) and Young et al. (2001). An overall SS model can be defined by these two equations: where x t is known as the state vector. In our case, the observation equation is defined in (5) and the state equation in (4). Then, x t = α 1,t , α * 1,t T . The remaining matrices and vectors are defined as follows: The white noise inputs, t , are assumed to be independent of the observation noise, σ 2 , and have a covariance matrix Q. These noise inputs determine the stochastic behaviour of the states. The SDP estimation can be divided into two steps. First, forward and backwards recursive algorithms are employed. Second, a backfitting algorithm is applied with the sorted data with respect to the dependent variable. Note that both the backwards pass smoothing and backfitting algorithms can be applied only if it is not necessary to work in real time.

Backwards pass smoothing equations
Note that Q r and P t are normalized with respect to the observation equation noise (σ 2 ), such as: P * t is the error covariance matrix related to the state estimates. The parameters inside that NVR matrix are estimated by a maximum likelihood prior to applying the recursive algorithms (Young et al. 2001).
Backfitting algorithm for SDP models.
1. Assume that FIS estimation has yielded an initial TVP estimate of α 0 1,t|N 2. Iterate k = 1, 2, … , k c a. sort both β A and β c according to the ascending order of β A . b. obtain an FIS estimate α k 1,t|N in the relationship a = 1 a ⋅ c (12)