key: cord-0887906-8gmjyzx8 authors: Rojas, Fernando; Wanke, Peter; Bravo, Fernando; Tan, Yong title: Inventory pooling decisions under demand scenarios in times of COVID-19 date: 2021-08-09 journal: Comput Ind Eng DOI: 10.1016/j.cie.2021.107591 sha: f98e5b8a357bdf582446c134e996fae5c647ca9b doc_id: 887906 cord_uid: 8gmjyzx8 Governments have been challenged to provide timely medical care to face the COVID-19 pandemic. The aim of this research is to propose a novel inventory pooling model to help determine order sizes and safety inventories in local hospital warehouses. The current study attempts to portray the availability of pharmaceutical items in public hospitals facing COVID-19 challenges. Different from previous studies, this research builds upon the consecrated theory of inventory pooling, extending it to pandemic circumstances where the intractability of kurtosis and skewness in inventory models are critical issues for making sure that medicines have high availability at a low cost. These effects on the total cost of inventory are explored and compared to a supply system with no consolidation. A continuous-review model is assumed with allocation rules for centralization and regular transshipment given different skewness and kurtosis structures for the demand, describing them by the copula method. This method models a multivariate demand considering that the marginal distributions of the demand can be specified by the Generalized Additive Model for Location, Scale and Shape, which offers advantages to model demands considering virtually any marginal statistical distribution. Numerical simulations and an illustrative example show that distributions of demands with more negative skewness and high kurtosis favor to a greater extent obtaining lower total costs with regular supply transshipment systems. Our study points out important considerations for supply chain decision makers when having demands with skewness and kurtosis patterns. The health care industry provides a vital service for modern societies (Kochan, Nowicki, Sauser, & Randall, 2018) . Analogously to what happens in other sectors, the health care industry competes based on time and quality. With the outbreak of COVID-19, several nations deployed a war effort in delivering goods and services to help control the pandemic, including the hospital care industry. In Chile, for instance, efforts have been made to increase the installed capacity for diagnoses, beds for critical patient care, as well as hospital equipment and supplies, including medication. Therefore, responsive and timely health care operations play a significant contribution to socio-economic integrity, promoting patient survival in rural and urban areas. What Newhouse (1970) said: ''hospital services seem to be desirable in some ethical sense, which justifies the claim that consumers have a right to medical care'', applies nowadays more than ever. We thus posit that the importance of hospitals in times of disaster provides leverage to consumer rights with respect to medical care, which in turn presents an indirect significant impact on the entire economy in the support and though there may not be a quantity of inventory on hand at the original supply location; see Evers (1996) and Rayat, Musavi, and Bozorgi-Amiri (2017) . When both systems are compared, the inventory pooling often provides lower costs than when using the IS of supply. The consolidated total cost (TC c ) of an inventory pooling and the total cost of an IS of supply (TC i ) are composed of various costs related to safety stocks and average quantities of inventories stored, ordered, and distributed during a cycle; see Knofius, van der Heijden, and Zijm (2019) and Wanke and Saliby (2009) . The TC c is subject to changes according to the allocation rules that determine if it is more convenient to centralize the warehouses to supply the demand or regularly occupy a transshipment policy; see Wen, Choi, and Chung (2019) . The decision of which allocation rule is more convenient for inventory pooling is made after comparing the corresponding TC c . To the best of our knowledge, the consolidation effect has been mostly studied assuming independent and normally distributed DPUTs, but when satisfying each point of sale or delivery, DPUTs are random variables (RVs) that may show any shape (Sadeghi & Niaki, 2015) . For the statistical distributions of RV independently and identically distributed (IID), such as the DPUTs , a statistical moment is a particular calculable dimension of the shape of its probability density functions (PDFs). The zero-th moment is always 1, the 1-th moment is the mean, the 2-th central moment is the variance, the 3-th standardized moment is the skewness, while the 4-th standardized moment is the kurtosis; see Casella and Berger (2002) , Deng, Miao, Ma, Wei, and Feng (2020) and Lin, Sun, and Yu (2020) . As we will see later in the background of this paper, all these moments are related, and for this paper we will postulate that they can influence inventory consolidation decisions. Often an RV data set has a joint statistical distribution based on marginal distributions with different skewness and kurtosis, as noted in Chan, Lim, and McAleer (2005) and Escribano and Pfann (1998) . This structure of dependence of the DPUT can be described by the copula method through the parametric specification and association of the marginal statistical distributions that will make up a multivariate joint distribution; see Autchariyapanitkul, Chanaim, and Sriboonchitta (2014) . Although theoretical foundations of copulas are complex, its practical treatment is made simple using open source software such as R (Kopczewski, Sobolewski, & Miernik, 2018) . This approach is made simpler than others, such as simulation (Qian, Li, & Hu, 2017) . In this context, it is essential to have excellent goodness of fit to actual data for a theoretical description of the marginal statistical distributions (Alavifard, 2019; Zhi, Wang, & Xu, 2020) . Generalized Additive Model for Location, Scale and Shape (GAMLSS) is a semiparametric regression type model introduced by Stasinopoulos and Rigby (2007) that allows great versatility in modeling random variables, which can be used in describing multivariate copulas of joint statistical distributions (Rohmer & Gehl, 2020) . In the current COVID-19 contingency scenario, drug demand patterns have undergone major changes. This is materialized in the fact that some products suffer large increases in the quantities demanded (demand shock), while others suffer an abrupt drop in these quantities requested for the treatment of patients. These changes in demand patterns can be shown as an increase in kurtosis and skewness to the right, in the first case, and as a decrease in kurtosis and left skewness in the second case. In both cases, these demand patterns are critical in determining the optimal lots to order and total inventory costs. In this study we extend the possibility of treating the kurtosis and asymmetry of the demand for critical items in the availability of pharmaceutical products in public hospitals that face COVID-19 contingency issues. In this context, the main contribution of this paper is to help understand how the kurtosis and skewness of the demand for items affect the supply systems of pooling and independent inventory policies and their total costs. So, the general objective of this paper is to propose a new methodology for studying how the skewness and kurtosis of joint distributions of DPUTs affect the TC c of an inventory pooling by using different allocation rules. The specific objectives are twofold: (i) compare the total costs between systems with inventory pooling and independent supply under a joint statistical distribution composed of marginal distributions with different moment measures; and (ii) analyze the changes in the indicators related to the inventory total costs under different DPUT scenarios with different moment measures. The remainder of the paper is organized as follows. In Section 2 we review literature regarding Healthcare inventory management and Inventory pooling with skewness and kurtosis scenarios in uncertain product demand. In the methodology of Section 3 we show how the multivariate statistical distribution of DPUTs based on copulas with different skewness and kurtosis of the marginal distributions affect the TC c of inventories with different allocation rules. Also in this section we design two simulation studies and an illustrative case with their respective statistical analyses. Results are shown in Section 4. Section 5 we provide managerial guidelines regarding the topic under study. Section 6 shows the discussion, limitations, and possible future research, and finally we presented our conclusions in Section 7. The literature has been investigating the issue of inventory management in the healthcare industry from different perspectives during recent years. Table 1 summarizes the related studies over the past decade. TC c usually corresponds to the sum of four components: (C1) consolidated safety stock (SS c ) multiplied by the holding cost (HC), (C2) consolidated cycle stock (CS c ) multiplied by HC, (C3) consolidated distribution cost (DC c ), and (C4) consolidated order cost (OC c ); see Wanke (2009) for details about these components. The uncertainty in the DPUT and lead time (LT) of the items is important in determining the sum of safety stocks that are consolidated in inventory pooling; see Tallon (1993) . Evers and Beier (1993) extended the SS c to multiple stocking points (locations or facilities) that serve the demand, while Wanke (2014) did the same with the CS c . In the context of reducing the TC c , Wanke (2009 Wanke ( , 2014 and Wanke and Saliby (2009) carried out a mathematical treatment of the consolidation effect considering three assumptions on the same safety-stock for the desired service level in all locations, which are as follows: 1. continuously review the model for lot-sizing; 2. define a reorder point as inventory control under uncertainty of the lead time demand (LTD), so that the consolidation does not affect the total average demand of the system; 3. variables DPUT, LT, and LTD must be independent and identically distributed in Gaussian form. The allocation rules for the demand are an important aspect to be considered in the consolidation effect. A first allocation rule corresponds to a given centralized location, which supplies the same proportion of demand to each decentralized location, known as the Tyagi and Das' allocation rule. Under this rule, and considering the DPUT as an IID RV, the maximum consolidation effect depends on the LT conditions between locations. Note that the consolidation effect is maximized when the TC c and inventory level are minimized. If both centralized locations present equal LT means and standard deviations (SDs), then the proportions of demand to be supplied in each decentralized location should be equal. Otherwise, that is under different LT means and/or SDs, a maximum consolidation effect is obtained when inventories are centralized into a single location; see Wanke (2009) . In a second allocation rule, a location can supply a primary demand area and also other demand areas, which in turn can be supplied by other locations. This is known as Ballou and Burnetas' allocation rule and here one must ask whether a proportion of the demand should Affiliation with local, regional, and national systems has mitigating effects under weak logistics services infrastructure with the mitigating effect being greatest for affiliation in local systems. Gebicki, Mooney, Chen, and Mazur (2014) Evaluate the medication inventory management Event-driven simulation The trade off between patient safety and cost can be addressed by incorporating drug characteristics in the ordering decisions. Wang, Cheng, Tseng, and Liu (2015) Examine hospital inventory management A dynamic drum-buffer-rope replenishment model The optimal replenishment timing and quantity of total inventory cost with no stock-out occurrence can be effectively determined by the model. Rahimi (2015) Examine the healthcare inventory routing problem A multi-objective mathematical model and possibilistic fuzzy approach The model proposed is superior in handling uncertain parameters. Forcina, Petrillo, Bona, Felice, and Silvestri (2017) Target stocking level which minimizes the total cost Dynamic models the system while satisfying the service level constraint The proposed model has characteristics of generality that allow the application in other areas Timajchi, Al-e Hashem, and Rekik (2019) Examine the healthcare inventory routing problem Bi-objective mixed integer mathematical programming Risk routes can be avoided and the economic supply network performance can be increased by the transshipment option. Saedi, Kundakcioglu, and Henry (2016) Investigate the optimal inventory policy for a healthcare facility A stochastic model Hospitals would benefit from inventory pooling since products in the healthcare industry, in particular in hospitals, have the characteristics that they have limited shelf-life. Chen, Xiao, Wang, and Lei (2020) Investigate the optimal order policies and the inventory adjustment quantity for a perishable products with a two-period half-life A stochastic optimization model 1. The expedited order plan performs better in terms of product wastage risk control, but it has a higher level of shortage risk compared to the returns plan. 2. The risk of wastage and shortage in inventory management can be effectively controlled by combining the expedited order plan and returns plan be supplied by a primary location and the remainder by secondary locations or not. In the assumption of IID DPUTs being considered, regular transshipment (RT) offers a good alternative for positively correlated DPUTs from all centralized locations with the possibility of balancing high/low values of the LT mean and DPUT SD at different centralized locations. It also seems to be the best system when the HC is very low; see Wanke and Saliby (2009) . As mentioned, previous works on the consolidation effect topic are based on the continuous review model under the assumption of IID normal DPUTs. These works provided evidence that the corresponding TC c decreased in different scenarios that considered the following indicators: (i) the correlation level between the DPUTs of the service points; (ii) the mean and SD of DPUTs and LTs; (iii) the security factor for the LTD; and (iv) holding, order, and distribution costs. For more details on these indicators and scenarios; see Wanke (2009) and Wanke and Saliby (2009) . We can point out the following main conclusions generated from previous works: (a) negative correlations between DPUTs and high values for SDs of DPUT are associated with a low TC c obtained by centralized systems; (b) positive correlations and low values for SDs of DPUT are associated with an IS of supply; (c) different SD (heterogeneity) for the LT is better handled by an RT because this system combines and balances the DPUTs that are served from different decentralized locations. The above findings are confirmed theoretically for non-Gaussian distributions by Corbett and Rajaram (2006) . However, our study differs from the previous one since we will explore the effect of left and right asymmetries and degree of kurtosis on DPUTs distributions to conclude how these influence the total costs between systems with inventory pooling and independent supply obtained by allocation rules. Skewness is a measure that reflects when a distribution or data set is symmetric or is shifted to the left or right of the central data. Kurtosis is a measure of data and outliers called tails, which can be heavy or light. Data with high kurtosis tend to have heavy tails or outliers, while a dataset with low kurtosis tends to have light tails or no outliers, where distributions with positive kurtosis are called leptokurtic. Those with kurtosis around zero are called mesokurtics and those with negative kurtosis are denominated platykurtics. Fig. 1 illustrates these differences. Note that the leptokurtic distributions concentrate higher probabilities of occurrence in the values close to the mean, while in the platykurtic the probabilities are more distributed in the tails of the distribution, while the mesokurtic distributions express the moderate cases of probabilistic distribution. On the other hand, the right skewness implies that there is a shift of the probabilities towards the data higher than the average, otherwise with the left skew, while the normal symmetry corresponds well to a Gaussian or normal distribution. Nelsen (1999) introduced the copula method to describe the dependence between RVs where the copula is a multivariate probability distribution with uniform marginal statistical distribution for each variable. Then, for the case of IID DPUTs, the copula method allows us to obtain the joint probabilistic distribution of correlated RVs without needing to know the joint distribution, which is easily obtainable from real DPUT data. As mentioned in the introduction, it is possible to use the GAMLSS model to model the marginal distributions of random variables that make up the joint multivariate distribution. Unlike generalized linear models (GLM), the GAMLSS considers a family of generalized, discrete, or continuous statistical distributions, which can have varying degrees of skewness and kurtosis. Thanks to its formulation, it is possible to model any parameter of these statistical distributions for the response variable linearly or not in an additive parametric or non-parametric form of covariates with known or random values; see Stasinopoulos, Rigby, and Akantziliotou (2008) . The advantage of using this type of statistical modeling is that GAMLSS is a regression toolbox appropriate for a big dataset of response variables that can consider linear or smoothing functions of predictive covariates to model any parameter of location, scale, or shape of the statistical distribution. The current packages available in R software (Stasinopoulos, Rigby, Voudouris, Heller, & De Bastiani, 2015) allow working with continuous (any type of skewness or kurtosis), discrete (including zero inflated data), and mixture statistical distributions. Models can be selected according to criteria of goodness of fit to the real data, as well as by generating random numbers with arbitrary distributions of interest for theoretical or empirical research (Rojas & Ibacache-Quiroga, 2020; Rojas, Leiva, Wanke, Lillo, & Pascual, 2019) . Our proposal differs from what has been addressed in the literature on the consolidation effect of inventories. We describe IID DPUT as a marginal statistical distribution described by GAMLSS models, making it possible to generate copulas to build joint multivariate with different skewness and kurtosis from the marginal distributions and to explore the TC c that leads to an inventory pooling by using different allocation rules. GAMLSS formulation. Let be an IID RV corresponding to the DPUT. If be the DPUT of an inventory item, we considered that is the expected value of a response variable. Consider to as a covariate. If ( | ) be a conditional PDF on parameters ( | is the conditional cumulative distribution function (CDF)), where = ( , , , ) ⊤ = ( 1 , 2 , 3 , 4 ) ⊤ is a vector of four distribution parameters. In the GAMLSS formulation, only is a function of the covariates, while and are location and scale parameters, and and are shape parameters. If { }, = 1, … , is an × 1 vector of the response variable to model, considering = 1, 2, 3, 4 as parameters, then is a link functions related to the th parameter and to covariates by the following additive models: where , , , , and 1 , for = 1, … , and = 1, 2, 3, 4, are × 1 vectors. 1 is an × 1 known matrix of variables and the regression coefficients 1 to be estimated is a 1 ×1 vector. ℎ is a semi-parametric additive function for the covariate evaluated at the vector , which is assumed fixed and known. For details of parameter estimate, diagnostic, and good fit on the data see Stasinopoulos and Rigby (2007) . Moments, skewness, and kurtosis. If is a CDF of any statistical distribution, considering the Riemann-Stieltjes integral (Liu, 2004) , the th moment of the statistical distribution is expressed by: with E as an expectation operator for the mean. The zero-th moment of any PDF is 1. The first raw moment is the mean: The second central moment is the variance, and the square root of the variance is the SD: The normalized th central moment of the RV is and represents the distribution. The normalized third central moment is called the skewness. A distribution that is skewed to the left has a negative skewness, and vice versa. Zero values indicate symmetry of the distribution. The Fisher coefficient of skewness (CSk) is defined as: where ′ 3 is the third moment centered. The fourth central moment is a measure of outliers values far from the average distribution values and is denominated kurtosis. Statistical distributions with kurtosis less than 3 are said to be ''platykurtic'', while distributions with kurtosis greater than 3 are said to be ''leptokurtic''. The Fisher coefficient of kurtosis (CK) is defined as: where ′ 4 is the fourth moment centered. Joint CDF 1 , 2 modeling by copulas. Now, let 1 and 2 be IID RVs with their unknown joint CDF 1 , 2 . Then, we have the relationship where is a parameter of dependence between 1 and 2 , which is denoted as = 1 , 2 . The copula defined in (5) is assumed to be continuous and twice differentiable. Thus, from (5) and Corbett and Rajaram (2006) , we obtain that the probability density function (PDF) associated with ( , , ), which can be expressed as and a quantile function (QF) ( ) = −1 ( ), for 0 < < 1. Note that −1 1 ( ) and −1 2 ( ), which can be obtained when marginal statistical distributions are known or adjusted by the GAMLSS model from actual data. Readers interested in using copulas to describe correlated DPUTs in consolidated effects are referred to Wanke (2014) . As mentioned in the introduction section, the TC c usually considers the components C1, C2, C3, and C4 that contain, respectively, the following elements: where, for the centralized location , OC is the order cost (expressed in USD$/order); HC = HC is the holding cost (expressed in USD$/unit); , is the proportion of the mean DPUT transferred from the decentralized location to the centralized location ; E( ) is the LT mean; Var( ) is the LT variance; and DC T , is the unitary transportation distribution cost to move a single item from a centralized location to a decentralized location . Note that 0 ≤ , ≤ 1, for all = 1, … , and = 1, … , , whereas ∑ =1 , = 1, for all , with being the number of centralized locations and the number of decentralized locations, for 1 ≤ ≤ . With and defined above, the case = = 2 can be extended to any value of , , with ≤ . This is based on the previous works of Ballou and Burnetas (2003) , Tyagi and Das (1998) , Wanke (2009 Wanke ( , 2014 and Wanke and Saliby (2009) , which allows us to have a benchmark for comparing our results. When = = 2, Tyagi and Das' allocation rule provides 1, = 2, = , for 0 ≤ ≤ 1, and ∑ =1 = 1, with 2 = 1 − 1 . However, under Ballou and Burnetas' allocation rule, it is known that 1,1 = 2,2 = , and 1,2 = 2,1 = 1 − , where the subindex denotes the primary location, that is, the location assigned with highest proportion of demand. Observe that the optimal solution under Ballou and Burnetas' allocation rule does not only behave differently from the solution obtained with Tyagi and Das' allocation rule, but it also implies in different inventory pooling models; see Wanke (2009) . Note that, under Ballou and Burnetas' allocation rule, both demand points are supplied exclusively by a warehouse, under an optimal solution where is zero or one (IS of supply), but under Tyagi and Das' allocation rule, if 1 is zero or one, both decentralized locations share one single serving facility. Nevertheless, when the optimal values of are greater than zero and less than one, RT takes place under Ballou and Burnetas' allocation rule. The above implies that both demand points are supplied by all the warehouses. 0 < 1 < 1 produces the same pattern as that Tyagi and Das' allocation rule, but they do not correspond to optimal solutions; see Wanke (2014) . The use of these allocation rules results in a consolidated inventory cost, allowing the decisions of both supply policies to be compared and to opt for the one with the lowest total cost. The differences and how these allocation rules operate are illustrated inFig. 2. Next we describe how to calculate CS c , SS c , OC c , and DC c under both allocation rules. Then, under Tyagi and Das' allocation rule in the decentralized location , expressions defined in (6)-(9) acquire the functional forms where , for = 1, 2, are defined as in (1), and In addition, to calculate DC c , we assume (without loss of generality) that inventories are consolidated at location 1. Under Ballou and Burnetas' allocation rule, and considering in the decentralized location , (6), (7), (8), and (9) acquire the following functional forms, respectively. ) , Note that, in (10), one must suppose that the highest proportions of are supplied with less DC T 1,1 and DC T 2,2 . In contrast, for the decentralized location considered under an IS of supply, the safety stock, cycle stock, order cost, and distribution cost are respectively given by Following Tallon (1993) , it is possible to define a consolidation effect indicator, which corresponds to ''the percentage reduction in aggregate safety stock made possible by consolidation effect of inventory from multiple locations into one location''. This indicator of consolidation effect of inventory into location is defined as As shown the SD of the DPUTs is an important parameter for determining both SS and SS as defined in (6) and (11), respectively. In turn, in Section 3.1 we have reviewed that the SD is an important input for calculating CSk and CK. Therefore, it is expected that these indicators influence SS and SS , and in turn the total costs reached in supply systems with inventory pooling under the allocation rules of Ballou and Burnetas (2003) and Tyagi and Das (1998) or independent systems. We present details of the computational framework utilized. We implemented our proposal in a non-commercial software named R; see http://www.r-project.org. See Rojas, Leiva, Wanke, and Marchant (2015) , Rojas, Wanke, Coluccio, Vega-Vargas, and Huerta-Canepa (2020), Wanke, Ewbank, Leiva, and Rojas (2016) and Wanke and Leiva (2015) to visualize R applications in supply models . The joint distribution when facing IID DPUT is simulated by using the copula method, occupying and R package denominated copula; see Hofert, Kojadinovic, Maechler, and Yan (2014) . The copula package provides classes of commonly used Archimedean, elliptical, and extreme value copulas, as well as other copula families such Gaussian. This package also contains methods for computing values related to PDF and CDF, generating a random number, plotting tools, fitting of copula models and goodness-of-fit tests. The mathematical programming of inventory pooling models is performed by an R package named nloptr; see Johnson (2014) . Recall that is the number of centralized locations and the number of decentralized locations, but the results presented here consider = = 2, because the results of this situation can be extended to any value of , , with ≤ . The simulation study 1 of scenarios establish different: (i) inventory policies (supply and allocation), (ii) statistical models for the DPUT, and (iii) mathematical programming to minimize costs. We used three types of supplies for inventory policies (IS/RT/inventory centralization -IC-) and two allocation rules (Tyagi and Das/Ballou and Burnetas). The statistical modeling is based on IID DPUTs assuming a Weibull type 3 (WEI3) statistical distribution, see Stasinopoulos et al. (2020) . We can estimate the parameters of the WEI statistical distribution proposed for the modeling of 1 and 2 , by using the histDist command of the gamlss package. We can also calculate the variances as: We assumed several structures of skewness and kurtosis with bivariate WEI3 marginal statistical distributions for the DPUT. We used the following indicators with WEI3 statistical distribution to generate these 10,000 scenarios: Statistical parameters • 1 , 2 ∼ U(80, 120), Mathematical programming is performed to minimize the total cost of inventory by using an objective function under Tyagi and Das' allocation rule, that is, when considering an IC. If the IS of supply is utilized, the mathematical programming for minimizing the total cost of inventory uses the objective function TC = HC(CS + SS ) + DC + OC , = 1, 2. Thus, such as Wanke (2009), we can compare IC and IS of supply based on (14) and (15), respectively, but now assuming WEI3 statistical distributions for the DPUTs. We call this Case 1. Furthermore, in this simulation study, we compare IC, RT, and IS of supply, such as in Wanke and Saliby (2009) , employing the total cost of inventory expressed as under both Tyagi and Das' and Ballou and Burnetas' allocation rules. Note that (16) is considered as the objective function to be minimized by the corresponding mathematical programming and compared to the optimization of (15). This is Case 2. Recall that the TC c defined in (14) is formed by the components C1, C2, C3, and C4, while, the TC c defined in (16) only considers components C1, C2, and C3. Analogously for the TC defined in (15), where the cost components considered are the same. Then, for = = 2, we considered 10,000 different scenarios to minimize the corresponding total cost, by using mathematical programming for obtaining the optimum values of C1, C2, C3, and C4. Our simulation study 2 is based on Wanke (2009) , who showed that the behavior of the consolidation effect indicator given in (13) can be treated as a response variable described by a linear function according to its relationship with the ratio = Var( 2 )∕Var( 1 ). Then, we explored the relationships of the consolidated effect in respect to different parameters of inventory model, for segments of . We illustrate our methodology with a real demand data set of two identical products with bivariate demand obtained from a mix of inventories that will allow us to show the methodology proposed. We first performed a descriptive statistical analysis for the demand data. Considering a goodness-of-fit analysis of the statistical distribution WEI3, we postulated it to theoretically describe the data and then estimated parameters that theoretically define their marginal statistical distribution and correlation. Second, we used inventory management models of centralization under the rules of Ballou and Burnetas (2003) and Tyagi and Das (1998) and compared their results of TCs with IS of supply, considering parameters, indicators, and the diverse costs involucrate. In this illustration, we performed a statistical analysis of the demand for two identical products used in palliative COVID-19 treatments in anonymous Chilean public hospital pharmacies. Our proposal is to be able to optimize an inventory system for these services and reduce costs associated with the centralization or decentralization of decisions. The data to be analyzed corresponded to a monthly demand of salbutamol inhalers from two hospital pharmacies. The actual supply system of this product is of the type IS. In simulation study 1, groups of scenarios were formed according to the grouping policy (IC, RT or IS) that achieves the lowest TC by allocation rules. We analyzed these results through the non-parametric Kruskall-Wallis test. This test makes it possible to relate the median of the different parameters of the inventory model used, in each of the groups (IC, RT, or IS) and to verify if these show statistically significant differences. For statistical analysis of simulation study 2, we explored the relationships of the consolidated effect with respect to different parameters of inventory model, through B-Spline regression, as this approach provides a way to estimate the values of a response variable between fixed points called knots. A polynomial regression between nodes is calculated, where the splines are a series of chained polynomial segments that are joined into knots. The choice of the degree of freedom (minus one if there is an intercept) generates knots at suitable quantiles of the variable, which will ignore missing values. An analogous regression analysis considering normal statistical distribution for DPUTs can be consulted in Table 7 of Wanke (2009) . We apply descriptive statistics to the DPUTs data for the statistical analysis of our illustrative actual case, then we estimate the parameters of the WEI 3 statistical distribution of these variables, by using the histDist command of the gamlss package in R software. We check the goodness of fit by examining their quantile-quantile (QQ) plots, which allows us to verify differences between the probability distribution of a population from which a random sample has been drawn and a theoretical distribution such as WEI3 used for comparison. Its interpretation is simple, since the graph should follow a diagonal within confidence intervals. Finally, by means of the optimization methods indicated in the preceding subsections, we conclude the best method of grouping inventories for the actual case. Firstly, in the 10,000 scenarios proposed for both Cases 1 and 2 and their skewness and kurtosis frameworks for the DPUTs, we analyzed which supply system would provide the minimum total cost. The results show that the total costs described by (14) for Case 1 or (16) for Case 2 are smaller when generated by an IC. For Case 1, we get 7590 optimum scenarios under IC and 2410 with an IS of supply. The number of scenarios with IC, RT, and IS of supply for Case 2 are 6270, 1830, and 1900, respectively. Secondly, using the 10,000 scenarios with IS, IC, and RT, we carried out a statistical analysis that is composed of two parts: (i) to identify the more relevant component (C1-C4) of the total cost from a quantitative perspective; and (ii) to determine whether statistically significant differences are presented for the association of the statistical indicators CSk, CK, and costs with the supply systems in Cases 1 and 2. See Fig. 3 to examine the shape of the statistical distribution and define the most important component of each total cost. Note that in all cases the values of the distribution median lead to the conclusion that C1 is the component with the highest proportion in total costs. We computed the number of scenarios where C1 is the larger component of total costs for Cases 1 and 2. For Case 1 with DPUTs, C1 corresponds to the highest proportion of the total cost in 8950 and 1050 scenarios under IC and IS of supply, respectively. For Case 2 with DPUTs, C1 is the most relevant component in 4570, 5310, and 120 scenarios under modalities RT, IC, and IS of supply, respectively. We compared the value of cost associated with component C1 regarding the sum of costs related to the remaining components (C2, C3, and C4) using the 2 test for the difference of proportions and detected statistically significant differences at 1% in favor of C1 for both Cases 1 and 2 under IC, RT, and IS of supply. This allows us to conclude that C1 is the most relevant component of the total cost. We used the Kruskal-Wallis (KW) test to compare medians of indicators related to statistical parameters, symmetry coefficient , kurtosis coefficient, and costs associated with IC, RT, and IS of supply. This test is carried out when minimum total costs are reached in Cases 1 (IC and IS of supply) and 2 (IC, RT and IS of supply). Tables 2 and 3 present the results for Cases 1 and 2, respectively. Note that the medians of 1 , 2 , 1 , 2 , CSk 1 , CSk 2 , CK 1 , and CK 2 are associated with a decrease in the total cost under Case 2, as indicated in Table 3 . Positive coefficients of correlation and small demand variability are associated with the IS of supply, confirming the findings of Wanke and Saliby (2009) . It is also possible to see that minors CSk 1 , CSk 2 , CK 1 , and CK 2 , are associated with the IS of supply, indicating that the skewness and kurtosis of the statistical distribution of the data of DPUTs are factors to consider in the decision of the best supply system (centralized or non-centralized). As pointed out in Section 3.1, these parameters directly impact to the variances, correlations, skewness, and kurtosis of DPUTs , which are elements of SS c and SS . Table 4 provides the medians of a selection of parameters that favor IC, RT, and IS of supply, considering only C1 or the HC of the safety stock (C1 = HC × SS). We realize that only the highest kurtosis and correlations favor lower C1 with RT systems and the effect of the other indicators is significant, but not always with the same type of association in relation to the total costs under (16), which balances the variances, skewness, and kurtosis of the better DPUTs. Furthermore, in comparison to total costs under (16), 1 , 2 does not need to be so negative for inducing lower C1 in an IC. However, it must increase to highly positive values to favor lower C1 in the IS of supply and to positive values close to zero to induce lower C1 in the RT. Table 5 shows the estimates of the regression coefficients obtained for the conditional expectance of the response variable related to consolidated effect (CE), standard errors (SE), and -value of the corresponding -test, with knots at suitable quantiles of < 1. The adjusted R-squared obtained for the B spline regression model was 0.7762. Note that not all the indicators turned out to be statistically significant as predictors of CE, and some of them changed the direction of the relationship by changing the knots at suitable quantile of CE related with . sd 1 , sd 2 and 1 , 2 , are always positive or negative predictors of CE, respectively, for both the indicators and their B splines. However, CSk 1 , CSk 2 , CK 1 and CK 2 , change the direction of the relationship with respect to CE according to the indicators or its B splines occupied as predictors. Table 3 Medians and corresponding Kruskal-Wallis -value when comparing them for the model mentioned, indicator, and supply system, based on the total cost for Case 2. The bivariate data set from the illustrative actual case has a length of 36 months (see Section 3.4), and its descriptive statistical measures can be consulted in Table 6 . Table 7 shows estimated parameters for the modeling by WEI3 statistical distribution of 1 and 2 by GAMLSS. In Fig. 4 , we show an empirical quantile-quantile (QQ) plot to check the good standing of our proposed statistical distribution (WEI3) to model the DPUT. Since almost every values is within the confidence bands, we conclude that the statistical distribution proposed fits the data perfectly. The rest of the coefficients required for carrying out the nonlinear optimizations under the centralization rules of Ballou and Burnetas (2003) and Tyagi and Das (1998) can be found in Table 8 with the results of TCs with respect to IS of supply compared. Now we can compare the TC obtained under the centralization rules of Ballou and Burnetas (2003) and Tyagi and Das (1998) with the IS supply system. Under the Tyagi and Das (1998) allocation rule, 1 = 0.999 (indicating centralization in localization 1) obtained a TC = 77.04 USD$. Under the Ballou and Burnetas (2003) allocation rule, = 0.624 (it is indicating an RT system) and obtained a TC = 99.42 USD$, while the IS of supply obtained a TC = 87.53 USD$, concluding that the best supply policy in this actual case is that of centralization in localization 1. Next we discuss the implications of our results for decision making management in the field of logistics by first referring to the importance of considering skewness and kurtosis patterns that DPUT data can have when inventory consolidation is desirable. With the results obtained, we conclude that these skewness and kurtosis patterns have a direct influence on parameters related to the variance and correlation of DPUTs. These parameters directly impact the component C1 of the total costs of inventory, which is the most relevant component in all the methods of inventory pooling and IS of supply; see Fig. 1 . A simple way to detect skewness and kurtosis patterns of the DPUT is by using histograms. For an interpretation of these graphs, the interested reader is referred to Brown (1997) . A selection of the main implications is provided next: • • In general, DPUT distributions with a more negative skewness (CSk < 0) and more leptokurtic, that is to say more pointed and with thicker tails than normal (CK > 3), favor to a greater extent obtaining lower total costs with RT systems. The RT of supply is favored given that it is associated with an increase in the variances and correlations of DPUTs in these conditions of symmetry and kurtosis. • In general, the above conditions of skewness and kurtosis favor a lower C1 with RT systems. On the contrary, more positive skewness of DPUTs and less leptokurtic, meaning tails not so heavy (CK tending to 3), favor to a greater extent the obtaining of a lower C1 with IS systems. • The statistical indicators with the highest direct impact on the percentages of consolidated effect are ordered as: DC T 2,2 > bs (sd 1 )1 > bs(CK 1 )2 > bs(sd 2 )1 > bs(sd 1 )2 > bs(sd 2 )2 > bs( 1 ) 1 > bs( 2 )2 > 2 . • The statistical indicators with the highest inverse impact on the percentages of consolidated effect are ordered as: DC T 1,2 > 1 , 2 > bs(CSk 1 )2 > bs(CK 2 )2 > bs( 1 , 2 )1 > HC > 2 > bs(DC T 2,2 )2. Fig. 5 shows a decision-making scheme developed to facilitate managing the supply decision of centralized or decentralized inventories according to skewness and kurtosis of the demand. To determine which type of inventory to pick, both skewness (CSk) and kurtosis (CK) must be evaluated simultaneously. The blue arrow indicates the direction to take that favors IC, the red arrow indicates the direction to take toward an IS, and the yellow arrow indicates the direction to take that favors RT. Coronavirus is not only a disease that has a devastating effect on people's lives by attacking the respiratory system, but it has a very fast spread rate, which has posed a big challenge to the related medicine production. The production process does not only need to focus on a timely medicine availability to meet a growing volume of demand, but also on how to simultaneously control the cost, which is the issue faced by all the producers. Since the outbreak of coronavirus, efforts have been made to investigate the issues of resource allocation, medical cost as well as supply of medical related resources (Bartsch et al., 2020; Zhang, Yao, Wang, Long, & Fu, 2020) . The above challenge was not only faced by the medical producers, but most importantly the hospitals that deal with this virus by directly contacting and treating the patients. Hospitals play a key role in slowing down the spread of the virus and there was a dilemma that was faced by the hospitals which was the significant increase in the number of patient's but limited amounts of resources in terms of medical related staffs, equipment, and medicine. Keeping these issues in mind, the investigation of pooling inventory policies would be of particularly relevant, useful, and important for policy making purposes by the government, through which a better decision-making process can be facilitated during the Covid-19 crisis. The aim of this study is to provide solutions to the policy makers in terms of how to mitigate inventory stock-out of critical supplies and significantly improve their ability by applying a novel stochastic inventory optimization model. Although this type of model has been proposed by a few studies in different contexts recently (Jackson, Tolujevs, & Reggelin, 2018; Pirhooshyaran & Snyder, 2020) , we are the first piece of research that uses this model in the context of COVID-19. The development of these practices would be of very paramount not only for the government as well as health care ministries in Chile, but it is very useful to develop a dynamic communication mechanism between hospitals transnationally, so it is recommended that the health managers should be actively engaged in this inventory policy approach and relevant trainings should be provided to them in terms of further improving their technical abilities of using R and Phyton, so as to better understand the policies. The usefulness of the current study also lies in the fact that the topic area is relevant not only to academic research in the area of operations research, but also very practical in the real world in terms of management of decisions. In the literature, there are a number of studies that investigate the issues of inventory control. Demand per unit of time is one of the areas of interest. Cobb, Rumí, and Salmerón (2013) investigated the optimal policies in inventory management when demand per unit of time follows a log-normal distribution. In comparison, while an inventory model is formulated by Pando, Garcı, San-José, Sicilia, et al. (2012) in which the demand per unit of time is a concave potential function of the inventory level. The stochastic inventory model was also examined by Hayya, Harrison, and Chatfield (2009) in which the demand per unit of time and the lead time random variables that are distributed independently and identically. Our results suggest new decision-making patterns such as considering the skewness and kurtosis of data related to demands per unit of time since it has shown its effect on the variance, which interacts with the best allocation rule for inventory pooling in terms of the consolidated safety stock multiplied by the holding cost and by the total cost. More precisely, in the presence of negative skewness and high kurtosis for the demands per unit of time, there must be fewer safety stocks and less total costs under the Ballou and Burnetas (2003) allocation rule that leads to inventory pooling for an RT system, implying that a contrary case favors dedicated facilities for the supply in terms of total cost. For the cost component related to the consolidated safety stock multiplied by the holding cost, the regular transshipment had a better statistical performance than the independent system of supply. This is in accordance with Reyes and Meade (2006) . We propose the methodology, which was implemented using the R software, facilitated by relevant R code . The methodology adopted in the current study possesses the advantages of being able to undertake a simulation study based on data with Skewness and Kurtosis. Regarding future research, this work is expandable to more general models with heteroscedasticity of variance, such as the widely known variants of generalized lineal models, which are very flexible allowing linear and non-linear functional structures. In addition, in line with this work, the methodology to model a time series, such as the generalized autoregressive and moving average (GARMA) model can be explored and its effects over CE. Theoretically, GARMA has been discussed by relevant studies (Albarracin, Alencar, & Ho, 2019; Gomes, Morettin, Cordeiro, & Taddeo, 2018) , however little effort has been made to apply this model to the real world practice, in particular in the area of inventory management. Furthermore, multivariate time series may be also considered. It is important to consider these new probabilistic approaches in cost reduction studies by potential locations for consolidation facilities that combine shipments, thus improving the level of service. Some of these issues are being analyzed by the authors whose findings will be reported in future articles. Due to the global outbreak of coronavirus at the beginning of this year, there has been an unprecedented pressure and challenge faced by hospitals across all countries. This also posed important questions in inventory management in the healthcare industry. Previously the empirical research focused on the cost perspective and in particular how to balance the cost and availability in the inventory management (Moons, Waeyenbergh, & Pintelon, 2019; Saedi et al., 2016) , while nowadays and in the future, we think the focus of the research in the area will be oriented to investigating how to manage the inventory during unexpected events. In other words, the future research trend will focus more on cost savings, inventory availability, uninterrupted supply, as well as on dealing with unexpected/unpredicted change in demand. Our study can be regarded as a pioneer in this research perspective by investigating inventory pooling decisions by considering the skewness and kurtosis of data. This research developed a novel inventory pooling model to assist in determining order sizes and safety inventories in public hospital warehouses, which have been dramatically impacted by COVID-19 challenges as regards availability of distinct items in distinct warehousing locations. Besides, the model developed here contributes to the current body of literature on inventory management, not only by overcoming the intractability of kurtosis and skewness in classic inventory models, but also by devising a framework to assess inventory consolidation gains. As long as all the moments of the distribution of the demand during the lead-time are related, and thus impact on inventory consolidation decisions, the copula method -through the parametric specification and association of the marginal statistical distributions that make-up a multivariate joint distribution -was employed for the first time to address inventory pooling issues. This paper is also helpful to hospital inventory managers and practitioners, ascertaining, under more realistic assumptions, that medicines have high availability at a low cost in different hospital warehouses. Therefore, in summary, our paper makes significant contributions from both theoretical and practical perspectives. Numerical simulations and an illustrative example were considered under the assumption of a continuous-review model while testing for distinct inventory allocation rules as regards inventory centralization and regular transshipment. Overall results indicated that distributions of demands with more negative skewness and high kurtosis favor are keys for obtaining lower total costs with regular supply transshipment systems. Based on the results, we also provided relevant practical policy implications, which was clearly presented in a good detail in the discussion section. Future studies should address not only distinct business environments as regards inventory cost parameters but should also consider different modeling assumptions, such as the periodic-review and the order-up-to-level regime. As regards the limitations of the proposed model, they are mainly related to modeling multivariate demand considering that the marginal distributions of the warehouse demands can be specified by the Generalized Additive Model for Location, Scale and Shape. While such approach offers advantages to model demands considering virtually any marginal statistical distribution, it falls short in addressing autocorrelated demands at each warehousing location. This could also constitute a research topic for future studies. Modelling default dependence in automotive supply networks using vine-copula Generalized autoregressive and moving average models: multicollinearity, interpretation and a new modified model Multi-commodity warehouse location and distribution planning with inventory consideration Portfolio optimization of stock returns in high-dimensions: A copula-based approach Estimating and auditing aggregate inventory levels at multiple stocking points Planning multiple location inventories The potential health care costs and resource use associated with covid-19 in the united states: A simulation estimate of the direct medical costs and health care resource use associated with covid-19 infections in the united states Collaborative management of inventory in Australian hospital supply chains: practices and issues. Supply Chain Management Skewness and kurtosis. Shiken: JALT Testing & Evaluation SIG Statistical inference Modelling multivariate international tourism demand and volatility Inventory strategies for perishable products with two-period shelf-life and lost sales Inventory management with log-normal demand per unit time A generalization of the inventory pooling effect to nonnormal dependent demand. Manufacturing and Service Operations Management The shaping of inventory systems in health services: A stakeholder analysis Reliability analysis of chatter stability for milling process system with uncertainties based on neural network and fourth moment method Systematic literature review on city logistics: overview, classification and analysis Optimising integrated inventory policy for perishable items in a multi-stage supply chain Non-linear error correction, asymmetric adjustment and cointegration The impact of transshipments on safety stock requirements The portfolio effect and multiple consolidation points: a critical assessment of the square root law An innovative model to optimise inventory management: a case study in healthcare sector Evaluation of hospital medication inventory policies Transformed symmetric generalized autoregressive moving average models Technical efficiency and its influencing factors in malaysian hospital pharmacy services A solution for the intractable inventory model when both demand and lead time are stochastic copula: Multivariate dependence with copulas: R package version 0.999-9 The combination of discrete-event simulation and genetic algorithm for solving the stochastic multi-product inventory optimization problem Consolidating spare parts for asset maintenance with additive manufacturing Impact of cloud-based information sharing on hospital supply chain performance: A system dynamics framework Bundling or unbundling? Integrated simulation model of optimal pricing strategies Behavioural data-driven analysis with bayesian method for risk management of financial services Refinement of an inequality of grüss type for riemann-stieltjes integral Measuring the logistics performance of internal hospital supply chains-a literature study An introduction to copulas Toward a theory of nonprofit institutions: An economic model of a hospital A multi-objective healthcare inventory routing problem; a fuzzy possibilistic approach Maximizing profits in an inventory model with both demand rate and holding cost per unit time dependent on the stock level Simultaneous decision making for stochastic multi-echelon inventory optimization with deep neural networks as decision makers A copula-based hybrid estimation of distribution algorithm for m-machine reentrant permutation flow-shop scheduling problem Bi-objective reliable locationinventory-routing problem with partial backordering under disruption risks: A modified amosa approach Improving reverse supply chain operational performance: a transshipment application study for not-for-profit organizations Sensitivity analysis of Bayesian networks to parameters of the conditional probability model using a beta regression approach A forecast model for prevention of foodborne outbreaks of non-typhoidal salmonellosis Modeling lot-size with time-dependent demand based on stochastic programming and case study of drug supply in chile Optimization of contribution margins in food services by modeling independent component demand Managing slow-moving item: a zero-inflated truncated normal approach for modeling demand Two parameter tuned multi-objective evolutionary algorithms for a bi-objective vendor managed inventory model with trapezoidal fuzzy demand Mitigating the impact of drug shortages for a healthcare facility: An inventory management approach Generalized additive models for location, scale and shape (GAMLSS) Instructions on how to use the gamlss package in r Flexible regression and smoothing: The gamlss packages in r. GAMLSS for Statistical Modelling. GAMLSS for Statistical Modeling The impact of inventory centralization on aggregate safety stock: the variable supply lead time case Inventory routing problem for hazardous and deteriorating items in the presence of accident risk with transshipment option Extension of the square root law for safety stock to demands with unequal variances Demand-pull replenishment model for hospital inventory management: a dynamic buffer-adjustment approach Consolidation effects and inventory portfolios Consolidation effects: Assessing the impact of tail dependence on inventory pooling using copulas Inventory management for new products with triangularly distributed demand and lead-time Exploring the potential use of the birnbaum-saunders distribution in inventory management Consolidation effects: Whether and how inventories should be pooled Fashion retail supply chain management: A review of operational models Supply chain risk management and hospital inventory: Effects of system affiliation Wuhan and hubei covid-19 mortality analysis reveals the critical role of timely supply of medical resources Impawn rate optimisation in inventory financing: A canonical vine copula-based approach This research was carried out thanks to the funding of the Fondecyt initiation project code: 11190004, Chile.