Estimation of Jointly Normally Distributed Demand for Cross-Selling Items in Inventory Systems with Lost Sales

Abstract

Demand estimation is often confronted with incomplete information of censored demand because of lost sales. Many estimators have been proposed to deal with lost sales when estimating the parameters of demand distribution. This study introduces the cross-selling effect into estimations, where two items are cross-sold because of the positive externality in a newsvendor-type inventory system. We propose an approach to estimate the parameters of a jointly normally distributed demand for two cross-selling items based on an iterative framework considering lost sales. Computational results based on more than two million numerical examples show that our estimator achieves high precision. Compared with the point estimations without lost sales, all the relative errors of the estimations of demand expectation, standard deviation, and correlation coefficient are no larger than 2% on average if the sample size is no smaller than 800. In particular, for demand expectation, the error is smaller than 1% if the comprehensive censoring level is no larger than four standard deviations (implying a $\mathrm{2}\sigma$-level of safety stock for each item), even if the sample size decreases to 50. This implies that the demand estimator should be competent in modern inventory systems that are rich in data.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In newsvendor-type inventory models, it is frequently assumed that the demand distribution for items is known or can be estimated based on sales data [1, 2]. Estimating the demand directly based on sales data implies that the sales quantity must be exactly equal to the realistic demand amount. However, this is not the case in practice [3]. The order quantity of an item is finite; thus, the inventory is limited but demand is often uncertain. It is very difficult to meet an uncertain demand based on a limited inventory. As a result, stock-outs emerge. Studies show that more than 85% of customers choose to give up the products when the items are out of stock [4], and therefore, the sales amount is the censored information of demand. When sales data is considered, demand might be underestimated. To remedy this flaw, demand estimation methods for censored sales data have been proposed [5, 6].

Another factor that should be considered in inventory management is the presence of externalities, including positive and negative ones [7]. As a type of positive externalities, cross-selling has not been considered in demand estimation. Cross-selling means that “a customer who has purchased a particular product may also be willing to purchase a related product” [8], and this implies that items are possibly purchased together owing to their unknown interior associations [9].

This study considers the cross-selling effect between two items to establish a new method to estimate the jointly (bivariate) normally distributed demand for two cross-selling items based on demand information from sales data that is censored by inventory order quantities. In our study, the sales data of any one item are not only censored because of lost sales caused by its own stock-out, but also affected by another associated item because of cross-selling. The core task is to recover the original demand distributions from the incomplete data that is a mixture of censored demand and cross sales. The proposed estimation is entirely based on the sales data, which is easily implemented in a modern business scenario. To the best of our knowledge, this is the first time that the demand estimation for cross-selling items with unobservable lost sales is proposed. Our work will consolidate the inventory model with positive externalities, as in the models of Netessine and Zhang [7] and Zhang et al. [10], and it can be used to provide basic information for inventory decision making.

The remainder of the paper is organized as follows. Section 2 reviews the related work. Section 3 includes assumptions and notations. Section 4 presents the demand estimation method and provides estimators. Section 5 illustrates the method by a numerical example, and Section 6 conducts numerical experiments to evaluate the performance of the estimators based on numerical examples. The paper is concluded in Section 7.

2. Literature Review

Many studies have focused on single-variable demand estimation considering censored data based on parametric methods where the type of random distribution is given. For censored data of a normal demand, Fisher [11] developed estimators for the mean and standard deviation by maximizing the likelihood function, where two specifically predetermined tables are needed. Cohen [12] introduced new formulae to improve Fisher's method, based on which the special tables were excluded. Later, Gupta [13] investigated the problem and presented a better linear unbiased estimator. Nahmias [5] developed a parametric methodology for estimating the normal demand and validated the approach based on censored samples. The method was further extended to the case of a negative binomial distribution [6]. Berk et al. [14] estimated the demand with censored data of negative binomial, gamma, Poisson, and normal distributions by Bayesian updating, adopting an approximate posterior distribution that matches the first two moments of the exact posterior. Heeseab [15] provided a Bayesian method for demand estimation considering unobserved lost sales, which depends on the prior information to estimate the posterior distribution. Jain et al. [16] modeled the cumulative demand arrivals during a period as a stochastic process with an unknown parameter. The parameter is updated in a Bayesian fashion by the observation of sales and stock-out events, which is further simplified by introducing stock-out timing observation.

Without a given type of demand distribution, nonparametric demand estimation methods are adopted. Beutel and Minner [17] proposed demand forecasting by solving a linear programming (LP) problem with fully observed demand, which does not need to know the random distribution. The method was further extended to deal with the censored observations of demand by Sachs and Minner [18]. They both assume that the demand is linearly dependent on multiple exogenous variables. For the nonlinear case, an intelligent forecasting algorithm based on artificial neural network (ANN) and conventional regression has been adopted. The method was applied for forecasting gasoline demand [19].

The above estimators are focused on the case of a single stochastic variable and did not consider consumption externalities. For multiple variables, Dempster et al. [20] proposed the expectation-maximization (EM) algorithm for the case of incomplete data, which, in contrast to the previous approaches, is an iterative method. Dahiya and Korwar [21] provided a parameter estimation method for the likelihood equation in the bivariate normal case. They assumed that the missing data were from only one variable and no consumption externalities existed. Adamids and Loukas [22] considered the missing values and proposed bivariate Poisson estimations for a two-item inventory based on the EM method without externalities. Regarding the time-varying demand, multivariate auto regressive integrated moving average (ARIMA) models were used to forecast the demand for perishable goods [23], but consumption externalities were not considered either.

Consumption externality means that the demand/consumption for one item affects the demand/consumption for other items, and there are two types of externality in inventory models: negative and positive [7]. Besanko et al. [24] presented a logit demand estimation based on the consumer choice model of substitutable items with negative externalities, where the endogenous price is the Nash equilibrium of wholesale and retail in a market. Anupindi et al. [25] and Smith and Agrawal [26] developed demand estimators considering stock-outs and multi-item substitutions, assuming Poisson arrivals of demand. Kök and Fisher [27] presented an EM-based method for estimating the parameters of substitution behavior and demand for products, which was used for optimizing the assortment planning of a retailer. Conlon and Mortimer [28] estimated the multinomial distribution demand under incomplete product availability by maximizing a likelihood function based on an EM algorithm. They assumed that the customer behavior can be described based on the discrete choice model of substitutable items. For multi-store multi-product substitution, Wan et al. [29] built demand estimation models considering that consumers may substitute a product for the stocked out item in the same store or switching to a neighboring store, where the parameters were estimated by the Markov chain Monte Carlo algorithm in a Bayesian manner.

It is noticeable that the extant inventory models and demand estimations with externalities are mainly focused on substitutable items, where choice happens because of substitution, and the customer behavior can be modeled by the consumer choice theory. However, this is not the case for cross-selling scenarios, where positive externalities exist between complementary items and no choice happens even if some items are stocked out. For comparing our study with previous works, Table 1 classifies the above-reviewed related literature into three groups.

Table 1

Summary of demand-estimation related literature.

Factors considered by estimation method	Literature
(1) Single item with stock-outs	Fisher [11], Cohen [12], Gupta [13], Agrawal and Smith [6], Berk et al. [14], Azadeh et al. [19], Heeseab [15], Beutel and Minner [17], Sachs and Minner [18], Jain et al. [16]
(2) Multiple items without externalities	Dempster et al. [20], Dahiya and Korwar [21], Adamids and Loukas [22], Huber et al. [23]
(3) Multiple items with externalities of substitution	Anupindi et al. [25], Besanko et al. [24], Smith and Agrawal [26], Kök and Fisher [27], Conlon and Mortimer [28], Wan et al. [29]

The related literature in Table 1 mainly includes three streams in the prior research: (1) estimation of the demand distribution of a single item considering censored data because of stock-out; (2) estimation of the jointly distributed demand for multiple items without considering externalities; and (3) estimation of the demand considering negative externalities because of substitution. We can observe that the extant research did not simultaneously consider the unobservable lost sales caused by stock-out and the positive externalities of cross-selling. Our study proposes a new method to estimate the jointly normally distributed demand for cross-selling items (because of positive externalities) based on censored sales data with unobservable lost sales.

3. Assumptions and Notations

This work is focused on inventory systems that are controlled based on the single-stage newsvendor strategy for multiple periods, but each period is independent. It means that neither the surplus product nor backordering is carried over from the former period to the next one. We focus on the cross-selling effect between two items, and thus the two-item inventory system is considered. Assumption 1 summarizes the above considerations.

Assumption 1.

The two-item inventory system is operated based on the single-stage newsvendor policy for each product, and there are no carryovers (surplus or backordered items) between adjacent periods.

Consider a two-item inventory system, where item $i$ and item $j$ affect each other in sales. In other words, besides the original “naked” demands for the two items (denoted by $D_i$ and $D_j$, respectively), the two items can also be demanded accompanying the sales of the counterpart item. For example, if a customer buys one unit of item $i$, she/he will possibly buy a certain quantity of the item $j$ because of the cross-selling effect with item $i$, and vice versa. Let ${\mu }_i,{\sigma }_i$ and ${\mu }_j,{\sigma }_j$ denote the mean and standard deviation of the demand for item $i$ and item $j$, respectively, and let $\rho =\mathrm{cov}(D_i,D_j)/{\sigma }_j{\sigma }_i$ represent the correlation coefficient. It is noted that externalities do not imply statistical correlations between the stochastic distributions of demand.

It could be considered that the cross-selling demand caused by the original “naked” demand might lead to the secondary “new” additional demand, which will make the problem too complex to solve. However, in practice, vendors mainly pay close attention to the cross-sold items caused by the original “naked” demand. The secondary “new” additional part is trivial and therefore is overlooked. Following this consideration, the following Assumption 2 holds.

Assumption 2.

Cross-selling demand is just caused by the original “naked” demand, and the cross-selling demand will not lead to a secondary “new” additional demand.

To quantitatively measure the cross-selling effect, we adopt the concept of association rules where cross-selling coefficients can be determined by data mining based on a transaction database [9, 30]. The cross-selling coefficient is defined as the amount of the cross-selling demand of one item caused by one unit another item's sale in average and vice versa. For clarity, we introduce Assumption 3.

Assumption 3.

Cross-selling coefficients are known and can be determined by mining of association rules, based on which the demand for an item has a linearly proportional increase by cross-selling.

Suppose that the sale of one unit of item $i$ (of the “naked” demand) leads to the sale of $r_{ij}$ units of item $j$, and the sale of one unit of item $j$ (of the “naked” demand) leads to the sale of $r_{ji}$ units of item $i$, in average. From Assumption 3, we can see that if there are no stock-outs, the total demand for the two items is $X_i=D_i+r_{ji}\bullet D_j$ and $X_j=D_j+r_{ij}\bullet D_i$, respectively.

Moreover, regarding the ordering policy of the inventory system, the following Assumption 4 is relevant in practice.

Assumption 4.

For each item, the order quantity is known and is constant for all periods.

Let $Q_i$ and $Q_j$ represent the order quantities of items $i$ and $j$, respectively. Denote $x_{it},x_{jt}$ the corresponding observed demands for item $i$ and item $j$ in period $t$, which are the realized sales that can be derived from transaction records. Evidently, the realized sales are bounded by order quantities, so that we have $x_{it}=\min (Q_i,X_{it})\le Q_i$ and $x_{jt}=\min (Q_j,X_{jt})\le Q_j$.

Note that if $x_{it}=Q_i$ or $x_{jt}=Q_j$, it could be considered that some sales are lost and the recorded sales amount is not equal to the actual demand. In this case, the sampled sales data may be censored as in Table 2.

Table 2

Sampled sales data being censored.

Period $t$

1

2

3

4

5

6

7

8

9

10

$$\dots$$

$$t$$

$$\dots$$

$T$

Sales of item $i$

$$x_{i1}$$

$$x_{i2}$$

$$Q_i$$

$$Q_i$$

$$x_{i5}$$

$$x_{i6}$$

$$x_{i7}$$

$$x_{i8}$$

$$Q_i$$

$$x_{i10}$$

$$\dots$$

$$x_{it}$$

$$\dots$$

$$x_{iT}$$

Sales of item $j$

$$Q_j$$

$$x_{j2}$$

$$x_{j3}$$

$$Q_j$$

$$x_{j5}$$

$$x_{j6}$$

$$Q_j$$

$$Q_j$$

$$x_{j9}$$

$$Q_j$$

$$\dots$$

$$x_{jt}$$

$$\dots$$

$$x_{jT}$$

It is noticeable that even if the original “naked” demand is smaller than the order quantity, this does not mean that there are no lost sales. Because of cross-selling, the sale of one item is influenced by another item's sales and order quantity. Suppose that the realized “naked” demand of items $i$ and $j$ in period 3 (see Table 2) is ${\xi }_{i3}$ and ${\xi }_{j3}$, respectively. Even ${\xi }_{i3}\le Q_i$, it is possible that $x_{i3}={\xi }_{i3}+r_{ji}\bullet \min \left(Q_j,{\xi }_{j3}\right)>Q_i$, resulting in lost sales of item $i$ because of the additional demand caused by cross-selling. As mentioned above, the cross-selling coefficients $r_{ji}$ and $r_{ij}$ have been determined by association rule mining, so we need to estimate the distribution parameters: ${\mu }_i$, ${\mu }_j$, ${\sigma }_i$, ${\sigma }_j$, and $\rho$, of the original demand $D_i$ and $D_j$. We summarize all notations in Table 3.

Table 3

Parameters and notations.

Parameters/notations	Meaning
$$i,j$$	Indices of items
$$D_i,D_j$$	Unknown original “naked” demand for items $i$ and $j$
$${\xi }_{it},{\xi }_{jt}$$	The realized original “naked” demand of items $i$ and $j$
$$x_{it},x_{jt}$$	Observed demand, i.e., the realized quantity of sales of items $i$ and $j$
$${\mu }_i,{\mu }_j$$	Means of demand (demand expectations) for items $i$ and $j$ to be estimated
$${\sigma }_i,{\sigma }_j$$	Standard deviation of demand for items $i$ and $j$ to be estimated
$$\rho$$	Correlation coefficient of the original demand to be estimated
$Q_i$ , $Q_j$	Order quantities of items $i$ and $j$
$$r_{ij},r_{ji}$$	Cross-selling coefficients
$${\stackrel{-}{x}}_i,{\stackrel{-}{x}}_j$$	Means of the sampled sales data of items $i$ and $j$
$${\tilde{\mu }}_i,{\tilde{\mu }}_j$$	Estimators of ${\mu }_i$ and ${\mu }_j$
$${\tilde{\sigma }}_i,{\tilde{\sigma }}_j$$	Estimators of ${\sigma }_i$ and ${\sigma }_j$
$$\tilde{\rho }$$	Estimator of $\rho$
$$\mathbb{z}/z$$	Quantile of standard normal distribution
$$T$$	Total number of sampled data (sample size)

4. Estimation Methodology

4.1. Estimators for One Item with Another Demand Given

Sorting increasingly all the sampled sales data, we will have two ordered data series as$$\begin{array}{}\text{(1)}& \mathrm{I}\mathrm{t}\mathrm{e}\mathrm{m}i:\left(x_{i\mathrm{1}},x_{i\mathrm{2}},\dots ,x_{iu},x_{i,u+\mathrm{1}},\dots ,x_{iT}\right),\text{(2)}& \mathrm{I}\mathrm{t}\mathrm{e}\mathrm{m}j:\left(x_{j\mathrm{1}},x_{j\mathrm{2}},\dots ,x_{jv},x_{j,v+\mathrm{1}},\dots ,x_{jT}\right),\end{array}$$where $x_{iu}<Q_i;x_{i,u+\mathrm{1}}=x_{i,u+\mathrm{2}}=\dots =x_{iT}=Q_i$ and $x_{jv}<Q_j;x_{j,v+\mathrm{1}}=x_{j,v+\mathrm{2}}=\dots =x_{jT}=Q_j$. $u$ and $v$ are the total number of $x_{it}$ and $x_{jt}$ that meet $x_{it}<Q_i$ and $x_{jt}<Q_j$, respectively.

Define$$\begin{array}{}\text{(3)}& {\stackrel{-}{x}}_u=\frac{\mathrm{1}}{u}{\displaystyle \sum }_{\tau =\mathrm{1}}^u{x}_{i\tau },\text{(4)}& {\stackrel{-}{x}}_v=\frac{\mathrm{1}}{v}{\displaystyle \sum }_{\tau =\mathrm{1}}^v{x}_{j\tau },\text{(5)}& S_u^{\mathrm{2}}=\frac{\mathrm{1}}{u-1}{\displaystyle \sum }_{\tau =\mathrm{1}}^u{\left(x_{i\tau }-{\stackrel{-}{x}}_u\right)}^{\mathrm{2}},\text{(6)}& S_v^{\mathrm{2}}=\frac{\mathrm{1}}{v-1}{\displaystyle \sum }_{\tau =\mathrm{1}}^v{\left(x_{j\tau }-{\stackrel{-}{x}}_v\right)}^{\mathrm{2}},\text{(7)}& P_i=\frac{u}{T},\text{(8)}& P_j=\frac{v}{T},\end{array}$$where ${\stackrel{-}{x}}_u$ and ${\stackrel{-}{x}}_v$ represent the mean values of the observed demand that are smaller than the order quantities, and $S_u^{\mathrm{2}},S_v^{\mathrm{2}}$ and $P_i,P_j$ are the corresponding sample variances and proportions, respectively. Regarding the distribution estimation of the original demand, we have the following results:

Proposition 5.

Given demand expectation ${\mu }_j$ and standard deviation ${\sigma }_j$ of item $j$, and with a known correlation coefficient $\rho$, the value of ${\mu }_i$ and ${\sigma }_i$ of item $i$ can be estimated by Eq. (9) and Eq. (10), respectively.$$\begin{array}{}\text{(9)}& {\tilde{\sigma }}_i=\sqrt{\left({\rho }^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}+\frac{S_u^{\mathrm{2}}}{\mathrm{1}-z_i\varphi \left(z_i\right)/P_i-{\varphi \left(z_i\right)}^{\mathrm{2}}/P_i^{\mathrm{2}}}}-{\rho r}_{ji}{\sigma }_j,\text{(10)}& {\tilde{\mu }}_i={\stackrel{-}{x}}_u+\frac{\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_i^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}+\mathrm{2}{\rho r}_{ji}{\tilde{\sigma }}_i{\sigma }_j}}{P_i}-r_{ji}{\mu }_j,\end{array}$$where $\varphi (\bullet )$ and $\Phi (\bullet )$ are the density function and cumulative probability function of the standard normal distribution, respectively, and $z_i={\Phi }^{-\mathrm{1}}(P_i)$ is the quantile of normal distribution corresponding to probability $P_i$.

Proof.

See Appendix A

Proposition 6.

Given ${\mu }_i$, ${\sigma }_i$, and $\rho$, the value of ${\mu }_j$ and ${\sigma }_j$ can be estimated by Eq. (11) and Eq. (12), respectively. $$\begin{array}{}\text{(11)}& {\tilde{\sigma }}_j=\sqrt{\left({\rho }^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\sigma }_i^{\mathrm{2}}+\frac{S_v^{\mathrm{2}}}{\mathrm{1}-z_j\varphi \left(z_j\right)/P_j-{\varphi \left(z_j\right)}^{\mathrm{2}}/P_j^{\mathrm{2}}}}-{\rho r}_{ij}{\sigma }_i,\text{(12)}& {\tilde{\mu }}_j={\stackrel{-}{x}}_v+\frac{\varphi \left(z_j\right)\sqrt{{\tilde{\sigma }}_j^{\mathrm{2}}+r_{ij}^{\mathrm{2}}{\sigma }_i^{\mathrm{2}}+\mathrm{2}{\rho r}_{ij}{\sigma }_i{\tilde{\sigma }}_j}}{P_j}-r_{ij}{\mu }_i,\end{array}$$where $z_j={\Phi }^{-\mathrm{1}}(P_j)$.

Proof.

It is straightforward from Proposition 5, because of the symmetry of the two items.

Evidently, the second terms at the right-hand sides of Eq. (9) and Eq. (11), $-{\rho r}_{ji}{\sigma }_j$ and $-{\rho r}_{ij}{\sigma }_i$, are both negative values for calculating the standard deviations of items $i$ and $j$. Not only that, the first terms under the radical sign, $({\rho }^2-1)r_{ji}^2{\sigma }_j^{\mathrm{2}}$ and $({\rho }^2-1)r_{ij}^2{\sigma }_i^{\mathrm{2}}$, are also negative. This implies that if we do not consider the cross-selling effect (the negative terms will not emerge), the standard deviation of the original demand might be overestimated. Similarly, the third terms at the right-hand sides of Eq. (10) and Eq. (12) also show that the cross-selling might lead to overestimation of the demand expectation. Thus, for a cross-selling inventory system, the association between items must be considered to mitigate the overestimation of the original demand. In other words, the cross-selling will increase either the observable demand amount or the observable demand uncertainty.

Besides, if we do not consider lost sales, the deviation and the mean of the original demand will be estimated based on $S_u^{\mathrm{2}}$, $S_v^{\mathrm{2}}$, ${\stackrel{-}{x}}_u$, and ${\stackrel{-}{x}}_v$. In this case, even if there is no cross-selling, the estimator will inevitably give underestimated outputs because of incomplete information from the censored sales data. This implies that lost sales decrease both the quantity and uncertainty of the effective demand.

4.2. Estimators of the Correlation Coefficient

To estimate the correlation coefficient $\rho$, define set $\mathcal{q}=\{(x_{i\tau },x_{j\tau })|x_{i\tau }<Q_i\mathrm{a}\mathrm{n}\mathrm{d}{x}_{j\tau }<Q_j\}$ representing all the sample points whose sales amounts are smaller than the order quantities. Evidently, for $(x_{i\tau },x_{j\tau })\in \mathcal{q}$, the sales data are not censored so that we have$$\begin{array}{c}\text{(13)}& x_{i\tau }={\xi }_{i\tau }+{\lambda }_{ji}{\xi }_{j\tau }{x}_{j\tau }={\xi }_{j\tau }+{\lambda }_{ij}{\xi }_{i\tau },\end{array}$$where ${\xi }_{i\tau }$ and ${\xi }_{j\tau }$ are the realized “naked” demand for items $i$ and $j$, respectively. Solving the above bivariate linear equation group yields$$\begin{array}{c}\text{(14)}& {\xi }_{i\tau }=\frac{x_{i\tau }-{\lambda }_{ji}{x}_{j\tau }}{\mathrm{1}-{\lambda }_{ij}{\lambda }_{ji}}{\xi }_{j\tau }=\frac{x_{j\tau }-{\lambda }_{ij}{x}_{i\tau }}{\mathrm{1}-{\lambda }_{ij}{\lambda }_{ji}}.\end{array}$$As the correlation is independent of the position where sample points are distributed, we will use ${\{{\xi }_{i\tau },{\xi }_{j\tau }\}}_{\tau \in \mathcal{q}}$ to estimate the correlation coefficient. Let $\left|\mathcal{q}\right|$ be the cardinality of set $\mathcal{q}$, and it follows that the estimator of the correlation coefficient is$$\begin{array}{}\text{(15)}& \tilde{\rho }=\frac{\mathrm{1}}{\left|\mathcal{q}\right|-\mathrm{1}}{{\displaystyle \sum }}_{\tau \in \mathcal{q}}\left(\frac{{\xi }_{i\tau }-{\stackrel{-}{\xi }}_{i\mathcal{q}}}{S_{i\mathcal{q}}}\right)\left(\frac{{\xi }_{j\tau }-{\stackrel{-}{\xi }}_{j\mathcal{q}}}{S_{j\mathcal{q}}}\right),\end{array}$$where ${\stackrel{-}{\xi }}_{i\mathcal{q}}$ and ${\stackrel{-}{\xi }}_{j\mathcal{q}}$ are the means of ${\{{\xi }_{i\tau }\}}_{\tau \in \mathcal{q}}$ and ${\{{\xi }_{j\tau }\}}_{\tau \in \mathcal{q}}$, respectively, and $S_{i\mathcal{q}}=\sqrt{(\mathrm{1}/(|\mathcal{q}|-\mathrm{1})){\sum }_{\tau \in \mathcal{q}}{\left({\xi }_{i\tau }-{\stackrel{-}{\xi }}_{i\mathcal{q}}\right)}^{\mathrm{2}}}$ and $S_{j\mathcal{q}}=\sqrt{(\mathrm{1}/(|\mathcal{q}|-\mathrm{1})){\sum }_{\tau \in \mathcal{q}}{\left({\xi }_{j\tau }-{\stackrel{-}{\xi }}_{j\mathcal{q}}\right)}^{\mathrm{2}}}$ are the corresponding sample standard deviations.

4.3. Iterative Estimation for Two Items

Propositions 5 and 6 give the demand estimation of one of the two cross-selling items when the other item's demand and correlation coefficient are given. However, the other item's original demand is unknown before accomplishing the estimation, and thus a circular deadlock is created. Here, we present an iterative procedure to obtain the estimation of the demand.

Specifically, we first estimate the correlation coefficient $\rho$ by Eq. (15) and obtain the estimator $\tilde{\rho }$. Then, we estimate the two items' demand without considering cross-selling based on the sampled data. This step is easy to realize based on Eqs. (9)–(12) by setting $r_{ij}=r_{ji}=0$, which yields $$\begin{array}{c}\text{(16)}& {\tilde{\sigma }}_{i\mathrm{0}}=\frac{P_i{S}_u}{\sqrt{P_i^{\mathrm{2}}-P_i{z}_i\varphi \left(z_i\right)-{\varphi \left(z_i\right)}^{\mathrm{2}}}},{\tilde{\mu }}_{i\mathrm{0}}={\stackrel{-}{x}}_u+\frac{\varphi \left(z_i\right){\tilde{\sigma }}_i}{P_i},\\ \text{(17)}& {\tilde{\sigma }}_{j\mathrm{0}}=\frac{P_j{S}_v}{\sqrt{P_j^{\mathrm{2}}-P_j{z}_j\varphi \left(z_j\right)-{\varphi \left(z_j\right)}^{\mathrm{2}}}},{\tilde{\mu }}_{j\mathrm{0}}={\stackrel{-}{x}}_v+\frac{\varphi \left(z_j\right){\tilde{\sigma }}_j}{P_j}.\end{array}$$Subsequently, substituting $({\tilde{\mu }}_{i0},{\tilde{\mu }}_{j0},{\tilde{\sigma }}_{i0},{\tilde{\sigma }}_{j0})$ and $\tilde{\rho }$ into Eqs. (9)–(12), we obtain the first iterative value of the estimators as$$\begin{array}{}\text{(18)}& {\tilde{\sigma }}_{i\mathrm{1}}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j\mathrm{0}}^{\mathrm{2}}+\frac{S_u^{\mathrm{2}}}{\mathrm{1}-z_i\varphi \left(z_i\right)/P_i-{\varphi \left(z_i\right)}^{\mathrm{2}}/P_i^{\mathrm{2}}}}-{\tilde{\rho }r}_{ji}{\tilde{\sigma }}_{j\mathrm{0}},\text{(19)}& {\tilde{\mu }}_{i\mathrm{1}}={\stackrel{-}{x}}_u+\frac{\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_{i\mathrm{1}}^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j\mathrm{0}}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ji}{\tilde{\sigma }}_{i\mathrm{1}}{\tilde{\sigma }}_{j\mathrm{0}}}}{P_i}-r_{ji}{\tilde{\mu }}_{j\mathrm{0}},\text{(20)}& {\tilde{\sigma }}_{j\mathrm{1}}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{i\mathrm{1}}^{\mathrm{2}}+\frac{S_v^{\mathrm{2}}}{\mathrm{1}-z_j\varphi \left(z_j\right)/P_j-{\varphi \left(z_j\right)}^{\mathrm{2}}/P_j^{\mathrm{2}}}}-{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{i\mathrm{1}}.\text{(21)}& {\tilde{\mu }}_{j\mathrm{1}}={\stackrel{-}{x}}_v+\frac{\varphi \left(z_j\right)\sqrt{{\tilde{\sigma }}_{j\mathrm{1}}^{\mathrm{2}}+r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{i\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{i\mathrm{1}}{\tilde{\sigma }}_{j\mathrm{1}}}}{P_j}-r_{ij}{\tilde{\mu }}_{i\mathrm{1}}.\end{array}$$Again, we substitute $({\tilde{\mu }}_{i1},{\tilde{\mu }}_{j1},{\tilde{\sigma }}_{i1},{\tilde{\sigma }}_{j1})$ and $\tilde{\rho }$ into Eqs. (9)–(12), yielding $({\tilde{\mu }}_{i2},{\tilde{\mu }}_{j2},{\tilde{\sigma }}_{i2},{\tilde{\sigma }}_{j2})$. Generally, for step $k\ge 1$, we have$$\begin{array}{}\text{(22)}& {\tilde{\sigma }}_{ik}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j,k-\mathrm{1}}^{\mathrm{2}}+\frac{S_u^{\mathrm{2}}}{\mathrm{1}-z_i\varphi \left(z_i\right)/P_i-{\varphi \left(z_i\right)}^{\mathrm{2}}/P_i^{\mathrm{2}}}}-{\tilde{\rho }r}_{ji}{\tilde{\sigma }}_{j,k-\mathrm{1}},\text{(23)}& {\tilde{\mu }}_{ik}={\stackrel{-}{x}}_u+\frac{\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_{ik}^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j,k-\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ji}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{j,k-\mathrm{1}}}}{P_i}-r_{ji}{\tilde{\mu }}_{j,k-\mathrm{1}},\text{(24)}& {\tilde{\sigma }}_{jk}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+\frac{S_v^{\mathrm{2}}}{\mathrm{1}-z_j\varphi \left(z_j\right)/P_j-{\varphi \left(z_j\right)}^{\mathrm{2}}/P_j^{\mathrm{2}}}}-{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{ik},\text{(25)}& {\tilde{\mu }}_{jk}={\stackrel{-}{x}}_v+\frac{\varphi \left(z_j\right)\sqrt{{\tilde{\sigma }}_{jk}^{\mathrm{2}}+r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{jk}}}{P_j}-r_{ij}{\tilde{\mu }}_{ik}.\end{array}$$The above procedures are repeated until a certain converging criterion is satisfied. We summarize the above procedure as follows.

Step 1.

Calculate $\tilde{\rho }$, the estimated value of correlation coefficient from Eq. (15).

Step 2.

Determine initial values of demand expectation and standard deviation for items $i,j$ from Eqs. (16) and (17).

Step 3.

For the $k$-th iteration ($k\ge 1$), estimate the demand expectation and standard deviation for items $i$ and $j$ from Eqs. (22)–(25).

Step 4.

When the relative changes of ${\mu }_i$, ${\mu }_j$, ${\sigma }_i$, and ${\sigma }_j$ between two adjacent iterations are smaller than $\delta$ for $q$ consecutive iterations, stop; otherwise, repeat Step 3.

$\delta$ is a small positive real number and $q$ is a positive integer. This study sets $\delta ={10}^{-6}$ and $q=10$. Regarding the convergence of the above procedure, we have Propositions 7 and 8.

Proposition 7.

For positive correlated demand, i.e., $\rho \ge 0$, the sequence $\left\{{\tilde{\sigma }}_{ik},{\tilde{\sigma }}_{jk}\right\}$ defined by Eqs. (22) and (24) are convergent.

Proof.

See Appendix B.

Proposition 8.

If the sequence $\left\{{\tilde{\sigma }}_{ik},{\tilde{\sigma }}_{jk}\right\}$ are convergent, the sequence $\left\{{\tilde{\mu }}_{ik},{\tilde{\mu }}_{jk}\right\}$ defined by Eqs. (23) and (25) is convergent.

Proof.

See Appendix B.

It is very difficult to theoretically prove the convergence in case of negative correlated demand $(\rho <0)$ because of the mathematical intractability. But the numerical experiment in Section 6 shows that, for each of more than two million numerical examples, the iterative procedures can quickly converge to a stable output even $\rho <0$. We conclude that, from a computational point of view, the iterative procedures based on our estimators are stable and convergent.

5. An Example for Illustration

5.1. The Sample Data

Suppose that the original “naked” demand for two cross-selling items is jointly distributed as $(D_i,D_j)~\mathcal{N}({\mu }_i,{\sigma }_i^2;{\mu }_j,{\sigma }_j^2;\rho )$, where ${\mu }_i=\mathrm{300},{\sigma }_i=\mathrm{10};{\mu }_j=\mathrm{400},{\sigma }_j=\mathrm{20};\rho =\mathrm{0.6}$. The two items are cross-sold with cross-selling coefficients $r_{ij}=0.2$ and $r_{ji}=0.3$. The order quantities are $Q_i=\mathrm{430},Q_j=\mathrm{470}$ in all periods.

If the original “naked” demand is realized as summarized in Table 4, the total demand without lost sales for the two items should be as summarized in Table 5, considering that $X_i=D_i+r_{ji}\bullet D_j$ and $X_j=D_j+r_{ij}\bullet D_i$.

Table 4

Original realized “naked” demand.

t	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
$${\xi }_{it}$$	313	311	289	303	292	296	290	272	293	286	305	294	308	315	283
$${\xi }_{jt}$$	395	418	391	411	420	377	388	363	369	356	412	422	412	425	367
t	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30
$${\xi }_{it}$$	309	290	299	299	299	291	295	284	286	292	301	309	300	295	303
$${\xi }_{jt}$$	404	366	390	396	375	418	397	388	375	418	382	422	392	386	387
t	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45
$${\xi }_{it}$$	299	307	317	290	301	289	306	306	282	305	301	287	286	290	321
$${\xi }_{jt}$$	387	428	416	378	405	410	437	402	374	401	406	392	406	380	421
t	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60
$${\xi }_{it}$$	294	269	285	314	302	306	305	309	295	288	321	321	305	291	297
$${\xi }_{jt}$$	381	350	380	419	420	403	396	382	395	386	401	438	428	369	393

Table 5

Total demand (without lost sales) for two items.

t	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
$$x_{it}$$	431.5	436.4	406.3	426.3	418.0	409.1	406.4	380.9	403.7	392.8	428.6	420.6	431.6	442.5	393.1
$$x_{jt}$$	457.6	480.2	448.8	471.6	478.4	436.2	446.0	417.4	427.6	413.2	473	480.8	473.6	488.0	423.6
t	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30
$$x_{it}$$	430.2	399.8	416.0	417.8	411.5	416.4	414.1	400.4	398.5	417.4	415.6	435.6	417.6	410.8	419.1
$$x_{jt}$$	465.8	424.0	449.8	455.8	434.8	476.2	456.0	444.8	432.2	476.4	442.2	483.8	452.0	445.0	447.6
t	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45
$$x_{it}$$	415.1	435.4	441.8	403.4	422.5	412.0	437.1	426.6	394.2	425.3	422.8	404.6	407.8	404.0	447.3
$$x_{jt}$$	446.8	489.4	479.4	436.0	465.2	467.8	498.2	463.2	430.4	462.0	466.2	449.4	463.2	438.0	485.2
t	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60
$$x_{it}$$	408.3	374.0	399.0	439.7	428.0	426.9	423.8	423.6	413.5	403.8	441.3	452.4	433.4	401.7	414.9
$$x_{jt}$$	439.8	403.8	437.0	481.8	480.4	464.2	457.0	443.8	454.0	443.6	465.2	502.2	489.0	427.2	452.4

Because of the limitation of order quantities, the observed demand (sales amount) will be censored as summarized in Table 6 according to $x_{it}=\min (Q_i,X_{it})$ and $x_{jt}=\min (Q_j,X_{jt})$.

Table 6

Observed demand censored by order quantities.

t	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
$$x_{it}$$	430	430	406.3	426.3	418.0	409.1	406.4	380.9	403.7	392.8	428.6	420.6	430	430	393.1
$$x_{jt}$$	457.6	470	448.8	470	470	436.2	446.0	417.4	427.6	413.2	470	470	470	470	423.6
t	16	17	18	19	20	21	22	23	24	25	26	27	28	29	30
$$x_{it}$$	430	399.8	416.0	417.8	411.5	416.4	414.1	400.4	398.5	417.4	415.6	406	430	419	424
$$x_{jt}$$	465.8	424.0	449.8	455.8	434.8	470	456.0	444.8	432.2	470	452	442.2	452.0	445.0	447.6
t	31	32	33	34	35	36	37	38	39	40	41	42	43	44	45
$$x_{it}$$	415.1	430	430	403.4	422.5	412.0	430	426.6	394.2	425.3	422.8	404.6	407.8	404.0	430
$$x_{jt}$$	446.8	470	470	436.0	465.2	467.8	470	463.2	430.4	462.0	466.2	449.4	463.2	438.0	470
t	46	47	48	49	50	51	52	53	54	55	56	57	58	59	60
$$x_{it}$$	408.3	374.0	399.0	430	428.0	426.9	423.8	423.6	413.5	403.8	430	430	430	401.7	414.9
$$x_{jt}$$	439.8	403.8	437.0	470	470	464.2	457.0	443.8	454.0	443.6	465.2	470	470	427.2	452.4

We will estimate ${\mu }_i,{\sigma }_i,{\mu }_j,{\sigma }_j$, and $\rho$ based on the observed sales amount in Table 6.

5.2. The Results

In total, there are 60 sampled data, among which 39 points meet $x_{it}<Q_i$ and $x_{jt}<Q_j$, and thus the coefficient $\rho$ is estimated from Eq. (15) as $\tilde{\rho }=0.5913$. Moreover, there are 14 data points of item $i$ censored and 18 data points of item $j$ censored so that $u=46$ and $v=42$. Thus, $P_i=u/T=46/60=0.7667$ and $P_j=v/T=42/60=0.7000$. We, therefore, have $z_i={\Phi }^{-\mathrm{1}}\left(P_i\right)={\Phi }^{-\mathrm{1}}\left(0.6833\right)=0.7279$ and $z_j={\Phi }^{-\mathrm{1}}\left(P_j\right)={\Phi }^{-\mathrm{1}}\left(0.6000\right)=0.5244$. Subsequently, we determine the initial value of $({\tilde{\mu }}_i,{\tilde{\sigma }}_i,{\tilde{\mu }}_j,{\tilde{\sigma }}_j)$ from Eqs. (16) and (17), yielding$$\begin{array}{c}\text{(26)}& {\tilde{\mu }}_{i0}=\mathrm{305.0}\mathrm{3},{\tilde{\sigma }}_{i0}=\mathrm{12.00};\\ {\tilde{\mu }}_{j0}=\mathrm{3}\mathrm{74.51},\\ {\tilde{\sigma }}_{j0}=\mathrm{20.88}.\end{array}$$Starting the iterative procedure, we will obtain the estimated results when the iteration converges, which yields$$\begin{array}{c}\text{(27)}& {\tilde{\mu }}_i=\mathrm{298.41},{\tilde{\sigma }}_i=\mathrm{1}\mathrm{2.01};\\ {\tilde{\mu }}_j=\mathrm{396.76},\\ {\tilde{\sigma }}_j=\mathrm{20.86}.\end{array}$$Figure 1 presents the converging process of the output by the estimators on each round of iteration.

[figure omitted; refer to PDF]

In Figure 1, the lines with triangles represent the point estimations obtained based on the sampling original “naked” demand $\left\{{\xi }_{it},{\xi }_{jt}\right\}$ in Table 4. It is easy to find that the point estimations are$$\begin{array}{c}\text{(28)}& {\tilde{\mu }}_{\xi i}=\frac{\mathrm{1}}{T}\sum {\xi }_{it}=\mathrm{298.18},{\tilde{\sigma }}_{\xi i}=\sqrt{\frac{\mathrm{1}}{T-\mathrm{1}}\sum {\left({\xi }_{it}-{\stackrel{-}{\xi }}_i\right)}^{\mathrm{2}}}=\mathrm{1}\mathrm{1.26},\\ {\tilde{\mu }}_{\xi j}=\frac{\mathrm{1}}{T}\sum {\xi }_{jt}=\mathrm{396.7}\mathrm{7},\\ {\tilde{\sigma }}_{\xi j}=\sqrt{\frac{\mathrm{1}}{T-\mathrm{1}}\sum {\left({\xi }_{jt}-{\stackrel{-}{\xi }}_j\right)}^{\mathrm{2}}}=\mathrm{20.41}.\end{array}$$Evidently, the results of our method are very close to the corresponding point estimations. The relative error is smaller than 1%. It shows that the performance of our estimators is as good as that of the point estimators in terms of the data sampled from the original demand. Similarly, we estimate $\rho$ based on $\left\{{\xi }_{it},{\xi }_{jt}\right\}$ and obtain ${\tilde{\rho }}_{\xi }=0.6557$. Our estimator gives the estimation $\tilde{\rho }=0.5913$, which is also close to ${\tilde{\rho }}_{\xi }$.

It is noticeable that we in fact cannot directly observe the realized original demand $\left\{{\xi }_{it},{\xi }_{jt}\right\}$ in practice. However, for the numerical experiments, we can use the point estimator based on the $\left\{{\xi }_{it},{\xi }_{jt}\right\}$ data as a benchmark to evaluate the accuracy of the new estimators proposed herein.

6. Numerical Experiments

We conduct numerical experiments with two aims: (1) to evaluate the estimation accuracy of the proposed method by comparing our estimators with the point estimation of the sampled original demand $\left\{{\xi }_{it},{\xi }_{jt}\right\}$; (2) to determine in a sensitivity analysis how the parameters impact the estimation.

6.1. Scheme for Generating Numerical Examples

The numerical experiment will be carried out based on artificially generated data. To investigate as many cases as possible, we will rearrange the values of the parameters from low to high, including the expectation, coefficient of variation (CV), and correlation coefficient of demand distribution, along with the cross-selling coefficient, censoring levels (order quantities), and sample size. Specifically, we select values of parameters for generating numerical examples, as summarized in Table 7.

Table 7

Selected values for parameters.

Parameter	Notation	Values to be selected
Expectation of item $i$’s original demand	$${\mu }_i$$	100, 300, 500
Expectation of item $j$’s original demand	$${\mu }_j$$	100, 300, 500
Coefficient of variation of item $i$	$${CV}_i={\sigma }_i/{\mu }_i$$	0.05, 0.15, 0.25
Coefficient of variation of item $j$	$${CV}_j={\sigma }_j/{\mu }_j$$	0.05, 0.15, 0.25
Cross-selling coefficient from item $i$ to $j$	$$r_{ij}$$	0.1, 0.3, 0.5, 0.7, 0.9
Cross-selling coefficient from item $j$ to $i$	$$r_{ji}$$	0.1, 0.3, 0.5, 0.7, 0.9
Correlation coefficient of demand	$$\rho$$	−0.9, −0.7, −0.5, −0.3, −0.1, 0.1, 0.3, 0.5, 0.7, 0.9
Order quantities (censoring level) of $i$	$$Q_i$$	${\mu }_{X_i}$ , ${\mu }_{X_i}+{\sigma }_{X_i}$, ${\mu }_{X_i}+\mathrm{2}{\sigma }_{X_i}$, ${\mu }_{X_i}+\mathrm{3}{\sigma }_{X_i}$
Order quantities (censoring level) of $j$	$$Q_j$$	${\mu }_{X_j}$ , ${\mu }_{X_j}+{\sigma }_{X_j}$, ${\mu }_{X_j}+\mathrm{2}{\sigma }_{X_j}$, ${\mu }_{X_j}+\mathrm{3}{\sigma }_{X_j}$
Sample size	$$T$$	50, 100, 300, 500, 700, 900, 1000

Note that ${\mu }_{X_i},{\mu }_{X_j}$ and ${\sigma }_{X_i},{\sigma }_{X_j}$ in Table 7 represent the expectation and standard deviation of total demand $X_i=D_i+r_{ji}\bullet D_j$ and $X_j=D_j+r_{ij}\bullet D_i$, respectively. According to Table 7, in total, we will test ${\mathrm{3}}^{\mathrm{2}}\times {\mathrm{3}}^{\mathrm{2}}\times {\mathrm{5}}^{\mathrm{2}}\times \mathrm{10}\times {\mathrm{4}}^{\mathrm{2}}\times \mathrm{7}=\text{2,268,000}$ numerical examples.

6.2. Computational Results

6.2.1. Relative Errors and Censoring Levels

As mentioned in Section 5.2, we will evaluate the performance of our estimation method by comparing the results with the point estimation of the data $\left\{{\xi }_{it},{\xi }_{jt}\right\}$ that is sampled from the original demand. We define the relative error of the estimated expectation of demand as$$\begin{array}{c}\text{(29)}& e_{\mu i}=\frac{\left|{\tilde{\mu }}_i-{\tilde{\mu }}_{\xi i}\right|}{{\tilde{\mu }}_{\xi i}}\times \mathrm{100}\%,e_{\mu j}=\frac{\left|{\tilde{\mu }}_j-{\tilde{\mu }}_{\xi j}\right|}{{\tilde{\mu }}_{\xi j}}\times \mathrm{100}\%,\end{array}$$where ${\tilde{\mu }}_{\xi i},{\tilde{\mu }}_{\xi j}$ are the corresponding point estimations. Similarly, the relative errors of the standard deviation and the correlation coefficient are$$\begin{array}{c}\text{(30)}& e_{\sigma i}=\frac{\left|{\tilde{\sigma }}_i-{\tilde{\sigma }}_{\xi i}\right|}{{\tilde{\sigma }}_{\xi i}}\times \mathrm{100}\%,e_{\sigma j}=\frac{\left|{\tilde{\sigma }}_j-{\tilde{\sigma }}_{\xi j}\right|}{{\tilde{\sigma }}_{\xi j}}\times \mathrm{100}\%,\\ e_{\rho }=\frac{\left|\tilde{\rho }-{\tilde{\rho }}_{\xi }\right|}{{\tilde{\rho }}_{\xi }}\times \mathrm{100}\%.\end{array}$$

Moreover, the censoring level is discriminated by the deviation of order quantities over the corresponding demand expectation, i.e., the multiples of ${\sigma }_{X_i}$ and ${\sigma }_{X_j}$. To be concise, we adopt the sum of the two multiples to comprehensively discriminate different censoring levels. The results will be reported considering the seven comprehensive censoring levels as those in the first column of Table 8.

Table 8

Comprehensive censoring levels.

Comprehensive censoring levels	Multiples of ${\sigma }_{X_i}$ and ${\sigma }_{X_j}$
$$\mathrm{0}\sigma$$	$$\left(0{\sigma }_{X_i},\mathrm{0}{\sigma }_{X_j}\right)$$
$$\mathrm{1}\sigma$$	$$\left(0{\sigma }_{X_i},1{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\left(1{\sigma }_{X_i},0{\sigma }_{X_j}\right)$$
$$\mathrm{2}\sigma$$	$$\left(0{\sigma }_{X_i},2{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\left(1{\sigma }_{X_i},1{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(2{\sigma }_{X_i},0{\sigma }_{X_j}\right)$$
$$\mathrm{3}\sigma$$	$$\left(0{\sigma }_{X_i},3{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\left(1{\sigma }_{X_i},2{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(2{\sigma }_{X_i},1{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(3{\sigma }_{X_i},0{\sigma }_{X_j}\right)$$
$$\mathrm{4}\sigma$$	$$\left(1{\sigma }_{X_i},3{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(2{\sigma }_{X_i},2{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(3{\sigma }_{X_i},1{\sigma }_{X_j}\right)$$
$$\mathrm{5}\sigma$$	$$\left(2{\sigma }_{X_i},3{\sigma }_{X_j}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(3{\sigma }_{X_i},2{\sigma }_{X_j}\right)$$
$$\mathrm{6}\sigma$$	$$\left(3{\sigma }_{X_i},3{\sigma }_{X_j}\right)$$

6.2.2. Results

For the seven censoring levels, we report the estimation errors of ${\mu }_i$, ${\mu }_j$, ${\sigma }_i$, and ${\sigma }_j$, under different sample sizes, correlation coefficients, cross-selling coefficients, and coefficients of variation, respectively.

(1) Errors vs. Sample Size ( $T$ ). Figure 2 presents the average error under different numbers of samples, i.e., $T=\mathrm{50,100,300,500,700,900,1000}$.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

We observe from Figure 2 that the demand expectation estimators yield very good estimation for the demand expectations of both items. Even for 50 sample points, the error is a maximum of 2.5%. With the standard deviation estimators, the error can reach more than 20% when the sample size is too small. However, if the sample size is larger than 800, the estimator performs much better, with an error no larger than 5%. We consider that for many modern business applications of inventory management, the data is no longer an obstruction. In particular, with the development of data acquisition and storage technologies, a large amount of data can be obtained from POS/MIS/MES, and big data systems [31]. The accuracy of our estimators can be further improved based on practical large sample data.

Moreover, it is expected that the estimation error will become smaller when censoring is alleviated. The results show that if the comprehensive censoring level reaches $\mathrm{4}\sigma$ (implying that the safety stock of the two items is in average $\mathrm{2}\sigma$-level with a fill-rate of 95.45%), the error of the demand expectation estimator decreases from about 2.0% of $\mathrm{0}\sigma$ to one no larger than 1.0%, and even to an error no larger than 0.5% with a $\mathrm{6}\sigma$ comprehensive censoring level (corresponding to a $\mathrm{3}\sigma$-level safety stock with a fill-rate of 99.74%). We also see that, with high fill-rate, the standard deviation error decreases from larger than 15% of small sample data to smaller than 2% of large sample data.

( 2) Errors vs. Correlation Coefficient ( $\rho$ ). For the value of demand expectation, the estimation error slightly increases with the correlation coefficient, except for the censoring level of $\mathrm{0}\sigma$. However, almost all estimation errors of the mean demand are no larger than 1.4%, and they even decrease to 0.2% when censoring is sufficiently alleviated with a comprehensive censoring level no smaller than $\mathrm{4}\sigma$ (corresponding to $\mathrm{2}\sigma$-level fill rate for each item in average), irrespective of the correlation coefficient.

It can also be observed that, under a higher safety stock level (with comprehensive censoring level no smaller than $\mathrm{4}\sigma$), the relative error of the standard deviation estimators is basically decreasing with the correlation coefficient. However, when censoring is aggravated (with censoring level no larger than $\mathrm{4}\sigma$), for weak correlation cases, the estimation error of the standard deviation reaches the maximal value, but it decreases when the correlation becomes stronger (irrespective of being positive or negative). The explanation is that, for a larger value of $\rho$, the estimator could capture more impact of the cross-selling effect (see Eqs. (9) and (11)), which might be beneficial to estimate the standard deviation of the demand.

( 3) Errors vs. Cross-Selling Coefficient ( $r_{ij},r_{ji}$ ). The cross-selling is of dual-direction, so that for each one of the seven censoring levels we must consider two impact factors: $r_{ij}$ and $r_{ji}$, which together will lead to a lengthy data report consisting of 28 two-dimensional tables. To be concise, we present a comprehensive index, $r=(r_{ij}+r_{ji})/\mathrm{2}$, to discriminate different levels of the cross-selling effect, as summarized in Table 9.

Table 9

Comprehensive cross-selling coefficients.

Comprehensive cross-selling coefficient $r=(r_{ij}+r_{ji})/\mathrm{2}$	Values of $(r_{ij},r_{ji})$
$$\mathrm{0.1}$$	$$(\mathrm{0.1,0.1})$$
$$\mathrm{0.2}$$	$$\left(\mathrm{0.1,0.3}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}(\mathrm{0.3,0.1})$$
$$\mathrm{0.3}$$	$$(\mathrm{0.1,0.5})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}(\mathrm{0.5,0.1})$$
$$\mathrm{0.4}$$	$$\left(\mathrm{0.1,0.7}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\left(\mathrm{0.3,0.5}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left(\mathrm{0.5,0.3}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.7,0.1})$$
$$\mathrm{0.5}$$	$$(\mathrm{0.1,0.9})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}(\mathrm{0.3,0.7})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.7,0.3})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.9,0.1})$$
$$\mathrm{0.6}$$	$$(\mathrm{0.3,0.9})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}(\mathrm{0.5,0.7})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.7,0.5})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.9,0.3})$$
$$\mathrm{0.7}$$	$$(\mathrm{0.5,0.9})\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.9,0.5})$$
$$\mathrm{0.8}$$	$$(\mathrm{0.7,0.9})\mathrm{o}\mathrm{r}(\mathrm{0.9,0.7})$$
$$\mathrm{0.9}$$	$$(\mathrm{0.9,0.9})$$

For the different comprehensive cross-selling coefficients in Table 9, the average relative errors of ${\mu }_i$, ${\mu }_j$, ${\sigma }_i$, and ${\sigma }_j$ are presented in Figure 4.

[figure omitted; refer to PDF]

The estimation errors of the demand expectation increase with the cross-selling coefficient, but they do not exceed 4.50% for all censoring levels. It means that the estimation for ${\mu }_i$ and ${\mu }_j$ has a wide range of applicability for handling different levels of cross-selling. For the standard deviation, the error tends to increase with the cross-selling coefficient, but it is still smaller than 5% for censoring levels of $\mathrm{5}\sigma$ and $\mathrm{6}\sigma$ in most cases. In particular, for $\mathrm{6}\sigma$-level censoring, the relative error is not really increasing with the cross-selling effect, and it even slightly decreases with larger cross-selling coefficients.

( 4) Errors vs. Coefficient of Variation $\left({CV}_i,{CV}_j\right)$ . The two coefficients of variation, ${CV}_i$ and ${CV}_j$, affect the estimation. Similarly, to be concise, we introduce a comprehensive coefficient of variation, $cv=({CV}_i+{CV}_j)/2$, to discriminate different variation levels, as summarized in Table 10.

Table 10

Comprehensive coefficient of variation.

Comprehensive coefficient of variation $cv=({CV}_i+{CV}_j)/2$	Values of $({CV}_i,{CV}_j)$
$$\mathrm{0.05}$$	$$(\mathrm{0.05,0.05})$$
$$\mathrm{0.1}$$	$$\left(\mathrm{0.05,0.15}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}(\mathrm{0.15,0.05})$$
$$\mathrm{0.15}$$	$$\left(\mathrm{0.05,0.25}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\left(\mathrm{0.15,0.15}\right)\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}\mathrm{\hspace{0.17em}\hspace{0.17em}}(\mathrm{0.25,0.05})$$
$$\mathrm{0.2}$$	$$(\mathrm{0.15,0.25})\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{o}\mathrm{r}(\mathrm{0.25,0.15})$$
$$\mathrm{0.25}$$	$$(\mathrm{0.25,0.25})$$

For the different comprehensive coefficients of variation in Table 10, the average relative errors of ${\mu }_i$, ${\mu }_j$, ${\sigma }_i$, and ${\sigma }_j$ are presented in Figure 5.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

Figure 5 presents that the estimation errors of ${\mu }_i$ and ${\mu }_j$ increase with the coefficient of variation, but even for a 20% variation both are still smaller than 1.5%. However, regarding the standard deviations, ${\sigma }_i$ and ${\sigma }_j$, the increase in order quantity for alleviating censoring is still an important means to obtain a good estimation. Only when the comprehensive censoring level reaches more than$\mathrm{5}\sigma$, the relative error becomes smaller than 5%.

(4) Estimation Errors of Correlation Coefficient ( $\rho$ ). We can see that the sample size is still a positive factor for the estimators of correlation coefficient, and it is evident that the more slightly the data is censored, the better the estimator will perform. We also observe that the error is smaller in negative correlations than in positive correlations on the overall trend. However, the estimation of $\rho$ still reaches a high precision provided that the fill rate of the inventory is high enough. Remarkably, for the $\mathrm{6}\sigma$ comprehensive censoring level, the error is no larger than 1.0% on average.

Combining all the above results, it is revealed that for inventory systems with at least $\mathrm{2}\sigma$-level safety stock (corresponding to $\mathrm{4}\sigma$ comprehensive censoring level), the estimation generates high-precision outcomes for larger sample data. It is worth noting that $\mathrm{2}\sigma$-level safety stock, i.e., 95% fill rate, is a basic service level for daily operations of a store. Hence, our approach could be of high practical relevance.

7. Conclusions

Demand distribution is key information for inventory decision making. Many estimators estimate the parameters of demand distribution for a single item based on the observed demand and considering unobservable lost sales. Methods focused on the demand estimation of multiple substitutable items have also been proposed. However, few studies consider the demand estimation for cross-selling items in inventory systems with lost sales, which leave a research gap in the extant literature.

This study extends the problem to the case of two cross-selling items. In the new scenario, the demand/sales of the two items are intertwined, and the lost sales may not be observed. We proposed a demand estimation method based on an iterative framework, which can deal with the cross-selling effect and lost sales. Numerical computations show that our estimators perform very well if the sample size is large enough. The methods will be competent in a modern inventory system that has substantial amounts of data and pursues a high-level service fill rate. The contribution of this study is summarized as follows:

(1) We propose an approach to estimate the jointly normally distributed demands for cross-selling items in inventory systems with lost sales. The approach can estimate the parameters of the original demand distribution from incomplete information that is a mixture of censored demand and cross sales.

(2) A closed-form estimator for an individual item's demand is developed, so that the parameters, including the means, deviations, and correlation coefficient of the cross-selling items' demands, can be efficiently calculated based on an iterative framework.

(3) The calculation of the proposed approach mainly utilizes the (censored) sales data, which should be easily implemented to provide demand information for inventory decisions in practical applications.

Our work can be further extended considering the following research lines. First, the method has not considered the case of multiple items. Second, we can further investigate different forms of demand distribution, for instance, Poisson distribution and negative binomial distribution. Third, the current method assumes that the expectation and standard deviation of the demand distribution are stable, but in many cases, the demand varies over time. Adapting the method to deal with a time-varying demand for cross-selling items will lead to more extensive researches.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The work is supported by National Natural Science Foundation of China (NSFC) under grant No. 71571006.

Appendix

Appendix

A.

Proof of Proposition 5.

( 1) Estimator of ${\mu }_i$ . Let $(D_i,D_j)~\mathcal{N}({\mu }_i,{\sigma }_i^2;{\mu }_j,{\sigma }_j^2;\rho )$ denote the original “naked” demand (jointly normally distributed) for items $i$ and $j$, where the correlation coefficient $\rho =\mathrm{cov}(D_i,D_j)/{\sigma }_j{\sigma }_i$, and we have$$\begin{array}{}\text{(A.1)}& D_i+r_{ji}{D}_j~\mathcal{N}\left({\mu }_i+r_{ji}{\mu }_j,{\sigma }_i^{\mathrm{2}}+\mathrm{2}{r}_{ji}\rho {\sigma }_i{\sigma }_j+r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}\right).\end{array}$$Let $f(x)$ denote the probability density function of a normally distributed variable $x~\mathcal{N}(\mu ,{\sigma }^2)$. Note that $E\left(D_i+r_{ji}{D}_j\right)={\mu }_i+r_{ji}{\mu }_j={\int }_{-\infty }^{+\infty }\left(D_i+r_{ji}{D}_j\right)f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)$, and it follows that$$\begin{array}{}\text{(A.2)}& E\left(D_i+r_{ji}{D}_j\right)={\mu }_i+r_{ji}{\mu }_j={{\displaystyle \int }}_{-\infty }^{Q_i}\left(D_i+r_{ji}{D}_j\right)f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)+{{\displaystyle \int }}_{Q_i}^{+\infty }\left(D_i+r_{ji}{D}_j\right)f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right).\end{array}$$Sort the sample points as Eqs. (1) and (2), and define$$\begin{array}{}\text{(A.3)}& {\stackrel{-}{x}}_u=\frac{\mathrm{1}}{u}{\displaystyle \sum }_{k=\mathrm{1}}^u\left(x_{ik}+r_{ji}{x}_{jk}\right),\end{array}$$where $x_{ik}$ and $x_{jk}$ are the recorded sales (the observed demand). Then, the first integral term in Eq. (A.2), i.e., the partial expectation, can be estimated by the sampled data as $$\begin{array}{}\text{(A.4)}& {{\displaystyle \int }}_{-\infty }^{Q_i}\left(D_i+r_{ji}{D}_j\right)f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)\approx \frac{\mathrm{1}}{T}{\displaystyle \sum }_{k=\mathrm{1}}^u\left(x_{ik}+r_{ji}{x}_{ik}\right)=\frac{u}{T}\left[\frac{\mathrm{1}}{u}{\displaystyle \sum }_{k=\mathrm{1}}^u\left(x_{ik}+r_{ji}{x}_{jk}\right)\right]=\frac{u}{T}{\stackrel{-}{x}}_u.\end{array}$$Defining $\mathbb{z}=(Q-\mu )/\sigma$, it is well known that$$\begin{array}{}\text{(A.5)}& {{\displaystyle \int }}_Q^{+\infty }xf\left(x\right)dx=\sigma \varphi \left(\mathbb{z}\right)+\mu \left[\mathrm{1}-\Phi \left(\mathbb{z}\right)\right],\end{array}$$where $\varphi$ and $\Phi$ are the probability density function and cumulative probability function of the standard normal distribution, respectively. Thus, for the second integral term in Eq. (A.2), we have$$\begin{array}{}\text{(A.6)}& {{\displaystyle \int }}_{Q_i}^{+\infty }\left(D_i+r_{ji}{D}_j\right)f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)=\sqrt{{\sigma }_i^{\mathrm{2}}+\mathrm{2}{r}_{ji}\rho {\sigma }_i{\sigma }_j+{r_{ji}^{\mathrm{2}}\sigma }_j^{\mathrm{2}}}\varphi \left(\mathbb{z}\right)+\left({\mu }_i+r_{ji}{\mu }_j\right)\left[\mathrm{1}-\Phi \left(\mathbb{z}\right)\right],\end{array}$$where $\mathbb{z}=(Q_i-\left({\mu }_i+r_{ji}{\mu }_j\right))/\sqrt{{\sigma }_i^{\mathrm{2}}+\mathrm{2}{r}_{ji}\rho {\sigma }_i{\sigma }_j+{r_{ji}^{\mathrm{2}}\sigma }_j^{\mathrm{2}}}$.

Substituting the above components into Eq. (A.2) yields$$\begin{array}{}\text{(A.7)}& \left({\mu }_i+r_{ji}{\mu }_j\right)\Phi \left(\mathbb{z}\right)=\frac{u}{T}{\stackrel{-}{x}}_u+\sqrt{{\sigma }_i^{\mathrm{2}}+\mathrm{2}{r}_{ji}\rho {\sigma }_j{\sigma }_i+r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}}\varphi \left(\mathbb{z}\right).\end{array}$$Note that $\Phi \left(\mathbb{z}\right)$ is the standard normal distribution probability that the stochastic variable is less than $Q_i$. Thus, we can estimate $\Phi \left(\mathbb{z}\right)$ by$$\begin{array}{}\text{(A.8)}& \Phi \left(\mathbb{z}\right)\approx P_i=\frac{u}{T},\end{array}$$and the estimator of $\mathbb{z}$ can be determined by $$\begin{array}{}\text{(A.9)}& \tilde{\mathbb{z}}={\Phi }^{-\mathrm{1}}\left(P_i\right)\end{array}$$Define$$\begin{array}{}\text{(A.10)}& z_i={\Phi }^{-\mathrm{1}}\left(P_i\right).\end{array}$$With the given estimator of ${\sigma }_i$ (denoted by ${\tilde{\sigma }}_i$), we estimate Eq. (A.7) as$$\begin{array}{}\text{(A.11)}& \left({\mu }_i+r_{ji}{\mu }_j\right)P_i=P_i{\stackrel{-}{x}}_u+\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_i^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}+\mathrm{2}\rho r_{ji}{\tilde{\sigma }}_i{\sigma }_j},\end{array}$$which yields the estimator of ${\mu }_i$ as $$\begin{array}{}\text{(A.12)}& {\tilde{\mu }}_i={\stackrel{-}{x}}_u+\frac{\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_i^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}+\mathrm{2}\rho r_{ji}{\tilde{\sigma }}_i{\sigma }_j}}{P_i}-r_{ji}{\mu }_j.\end{array}$$This is Eq. (10).

( 2) Estimator of ${\sigma }_i$ . Consider the variance of $D_i+r_{ji}{D}_j$, and we have $$\begin{array}{}\text{(A.13)}& {\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}={\sigma }_i^{\mathrm{2}}+\mathrm{2}{r}_{ji}\rho {\sigma }_i{\sigma }_j+r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}={{\displaystyle \int }}_{-\infty }^{Q_i}{\left(D_i+r_{ji}{D}_j-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)+{{\displaystyle \int }}_{Q_i}^{+\infty }{\left(D_i+r_{ji}{D}_j-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right).\end{array}$$Note that we can estimate the first integral term of Eq. (A.13) by$$\begin{array}{}\text{(A.14)}& S_T^{\mathrm{2}}=\frac{\mathrm{1}}{T-\mathrm{1}}{\displaystyle \sum }_{k=\mathrm{1}}^u{\left(x_{ik}+r_{ji}{x}_{jk}-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}.\end{array}$$Extending $S_T^2$ yields$$\begin{array}{}\text{(A.15)}& S_T^{\mathrm{2}}=\frac{\mathrm{1}}{T-\mathrm{1}}{\displaystyle \sum }_{k=\mathrm{1}}^u{\left(x_{ik}+r_{ji}{x}_{jk}-{\stackrel{-}{x}}_u+{\stackrel{-}{x}}_u-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}=\frac{\mathrm{1}}{T-\mathrm{1}}\left\{\begin{array}{c}{\displaystyle \sum }_{k=\mathrm{1}}^u\left[{\left(x_{ik}+r_{ji}{x}_{jk}-{\stackrel{-}{x}}_u\right)}^{\mathrm{2}}+\mathrm{2}\left(x_{ik}+r_{ji}{x}_{jk}-{\stackrel{-}{x}}_u\right)\left({\stackrel{-}{x}}_u-{\mu }_i-r_{ji}{\mu }_j\right)\right]\\ +u{\left({\stackrel{-}{x}}_u-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}\end{array}\right\}.\end{array}$$Because the middle term is zero, we obtain(A.16) S T 2 = 1 T - 1 ∑ k = 1 u x i k + r j i x j k - x - u 2 + u T - 1 x - u - μ i - r j i μ j 2 = u - 1 T - 1 S u 2 + u T - 1 x - u - μ i - r j i μ j 2 ,

where$$\begin{array}{}\text{(A.17)}& S_u^{\mathrm{2}}=\frac{\mathrm{1}}{u-\mathrm{1}}{\displaystyle \sum }_{k=\mathrm{1}}^u{\left(x_{ik}+r_{ji}{x}_{jk}-{\stackrel{-}{x}}_u\right)}^{\mathrm{2}}.\end{array}$$Thus, we can estimate the first integral term as$$\begin{array}{}\text{(A.18)}& {{\displaystyle \int }}_{-\infty }^{Q_i}{\left(D_i+r_{ji}{D}_j-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)=\frac{u-\mathrm{1}}{T-\mathrm{1}}{S}_u^{\mathrm{2}}+\frac{u}{T-\mathrm{1}}{\left({\stackrel{-}{x}}_u-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}.\end{array}$$Moreover, for a normally distributed variable $x~\mathcal{N}(\mu ,{\sigma }^2)$, it is well known that$$\begin{array}{}\text{(A.19)}& {{\displaystyle \int }}_Q^{+\infty }{\left(x-\mu \right)}^{\mathrm{2}}\varphi \left(x\right)dx={\sigma }^{\mathrm{2}}\left[\mathrm{1}-\Phi \left(\mathbb{z}\right)\right]+{\sigma }^{\mathrm{2}}\mathbb{z}\varphi \left(\mathbb{z}\right).\end{array}$$Therefore, for the second integral term of Eq. (A.13), we have$$\begin{array}{}\text{(A.20)}& {{\displaystyle \int }}_{Q_i}^{+\infty }{\left(D_i+r_{ji}{D}_j-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}f\left(D_i+r_{ji}{D}_j\right)d\left(D_i+r_{ji}{D}_j\right)={\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}\left[\mathrm{1}-\Phi \left(\mathbb{z}\right)\right]+{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}\mathbb{z}\varphi \left(\mathbb{z}\right).\end{array}$$Substituting Eqs. (A.18) and (A.20) into Eq. (A.13) yields$$\begin{array}{}\text{(A.21)}& {\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}=\frac{u-\mathrm{1}}{T-\mathrm{1}}{S}_u^{\mathrm{2}}+\frac{u}{T-\mathrm{1}}{\left({\stackrel{-}{x}}_u-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}+{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}\left[\mathrm{1}-\Phi \left(\mathbb{z}\right)\right]+{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}\mathbb{z}\varphi \left(\mathbb{z}\right),\end{array}$$i.e., $$\begin{array}{}\text{(A.22)}& \mathrm{0}=\frac{u-\mathrm{1}}{T-\mathrm{1}}{S}_u^{\mathrm{2}}+\frac{u}{T-\mathrm{1}}{\left({\stackrel{-}{x}}_u-{\mu }_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}-{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}\Phi \left(\mathbb{z}\right)+{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}\mathbb{z}\varphi \left(\mathbb{z}\right).\end{array}$$Similarly, we can estimate that $\Phi \left(\mathbb{z}\right)\approx P_i=u/T$. Then, based on the estimators ${\tilde{\mu }}_i$ for ${\mu }_i$, we reformulate Eq. (A.22) as$$\begin{array}{}\text{(A.23)}& \mathrm{0}=\frac{u-\mathrm{1}}{T-\mathrm{1}}{S}_u^{\mathrm{2}}+\frac{u}{T-\mathrm{1}}{\left({\stackrel{-}{x}}_u-{\tilde{\mu }}_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}-{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}{P}_i+{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}{z}_i\varphi \left(z_i\right)\end{array}$$where $z_i={\Phi }^{-\mathrm{1}}\left(P_i\right)$.

Moreover, according to Eq. (A.12), we have$$\begin{array}{}\text{(A.24)}& {\left({\stackrel{-}{x}}_u-{\tilde{\mu }}_i-r_{ji}{\mu }_j\right)}^{\mathrm{2}}=\frac{{\tilde{\sigma }}_{D_i+r_{ji}{D}_j}^{\mathrm{2}}{\varphi \left(z_i\right)}^{\mathrm{2}}}{{P_i}^{\mathrm{2}}}.\end{array}$$where ${\tilde{\sigma }}_{D_i+r_{ji}{D}_j}^{\mathrm{2}}={\tilde{\sigma }}_i^{\mathrm{2}}+2r_{ji}\rho {\tilde{\sigma }}_i{\sigma }_j+r_{ji}^2{\sigma }_j^2$ is estimated by the given estimator ${\tilde{\sigma }}_i$ of ${\sigma }_i$. Substituting Eq. (A.24) into Eq. (A.23) yields$$\begin{array}{}\text{(A.25)}& \mathrm{0}=\frac{u-\mathrm{1}}{T-\mathrm{1}}{S}_u^{\mathrm{2}}+\frac{u{\tilde{\sigma }}_{D_i+r_{ji}{D}_j}^{\mathrm{2}}{\varphi \left(z_i\right)}^{\mathrm{2}}}{\left(T-\mathrm{1}\right){P_i}^{\mathrm{2}}}-{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}{P}_i+{\sigma }_{D_i+r_{ji}{D}_j}^{\mathrm{2}}{z}_i\varphi \left(z_i\right),\end{array}$$and we obtain$$\begin{array}{}\text{(A.26)}& {\tilde{\sigma }}_{D_i+r_{ji}{D}_j}^{\mathrm{2}}=\frac{\left(\left(u-\mathrm{1}\right)/\left(T-\mathrm{1}\right)\right)S_u^{\mathrm{2}}}{P_i-z_i\varphi \left(z_i\right)-u{\varphi \left(z_i\right)}^{\mathrm{2}}/\left(T-\mathrm{1}\right)P_i^{\mathrm{2}}}.\end{array}$$Substituting ${\tilde{\sigma }}_{D_i+r_{ji}{D}_j}^{\mathrm{2}}={\tilde{\sigma }}_i^{\mathrm{2}}+2r_{ji}\rho {\tilde{\sigma }}_i{\sigma }_j+r_{ji}^2{\sigma }_j^{\mathrm{2}}$ into Eq. (A.26) and solving the quadratic equation with respect to ${\tilde{\sigma }}_i$ yields$$\begin{array}{}\text{(A.27)}& {\tilde{\sigma }}_i=\sqrt{\left({\rho }^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}+\frac{\left(u-\mathrm{1}\right)S_u^{\mathrm{2}}{P}_i^{\mathrm{2}}}{\left[P_i-z_i\varphi \left(z_i\right)\right]\left(T-\mathrm{1}\right)P_i^{\mathrm{2}}-u{\varphi \left(z_i\right)}^{\mathrm{2}}}}-r_{ji}\rho {\sigma }_j.\end{array}$$When $T$ and $u$ are large enough, $(u-\mathrm{1})/(T-\mathrm{1})$ and $u/(T-\mathrm{1})$ could be approximated by $P_i=u/T$; then, Eq. (A.27) is simplified as$$\begin{array}{}\text{(A.28)}& {\tilde{\sigma }}_i=\sqrt{\left({\rho }^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\sigma }_j^{\mathrm{2}}+\frac{S_u^{\mathrm{2}}}{\mathrm{1}-z_i\varphi \left(z_i\right)/P_i-{\varphi \left(z_i\right)}^{\mathrm{2}}/P_i^{\mathrm{2}}}}-\rho r_{ji}{\sigma }_j.\end{array}$$This is Eq. (9).

B.

Proof of Proposition 7.

We will prove that the sequences $\left\{{\tilde{\sigma }}_{ik}\right\}$ and $\left\{{\tilde{\sigma }}_{jk}\right\}$ are monotonic and bounded, so that they are convergent. To be concise, define $$\begin{array}{c}\text{(B.1)}& {\widehat{S}}_v^{\mathrm{2}}=\frac{S_v^{\mathrm{2}}}{\mathrm{1}-z_j\varphi \left(z_j\right)/P_j-{\varphi \left(z_j\right)}^{\mathrm{2}}/P_j^{\mathrm{2}}},{\widehat{S}}_u^{\mathrm{2}}=\frac{S_u^{\mathrm{2}}}{\mathrm{1}-z_i\varphi \left(z_i\right)/P_i-{\varphi \left(z_i\right)}^{\mathrm{2}}/P_i^{\mathrm{2}}}.\end{array}$$From Eqs. (22) and (24), we have$$\begin{array}{}\text{(B.2)}& {\tilde{\sigma }}_{j,k+\mathrm{1}}-{\tilde{\sigma }}_{jk}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{i,k+\mathrm{1}}^{\mathrm{2}}+{\widehat{S}}_v^{\mathrm{2}}}-\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+{\widehat{S}}_v^{\mathrm{2}}}-{\tilde{\rho }r}_{ij}\left({\tilde{\sigma }}_{i,k+\mathrm{1}}-{\tilde{\sigma }}_{ik}\right)\text{(B.3)}& {\tilde{\sigma }}_{i,k+\mathrm{1}}-{\tilde{\sigma }}_{ik}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{jk}^{\mathrm{2}}+{\widehat{S}}_u^{\mathrm{2}}}-\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j,k-\mathrm{1}}^{\mathrm{2}}+{\widehat{S}}_u^{\mathrm{2}}}-{\tilde{\rho }r}_{ji}\left({\tilde{\sigma }}_{jk}-{\tilde{\sigma }}_{j,k-\mathrm{1}}\right)\text{(B.4)}& {\tilde{\sigma }}_{jk}-{\tilde{\sigma }}_{j,k-\mathrm{1}}=\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+{\widehat{S}}_v^{\mathrm{2}}}-\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{i,k-\mathrm{1}}^{\mathrm{2}}+{\widehat{S}}_v^{\mathrm{2}}}-{\tilde{\rho }r}_{ij}\left({\tilde{\sigma }}_{ik}-{\tilde{\sigma }}_{i,k-\mathrm{1}}\right)\end{array}$$It is easy to obtain the convergence in case of ${\tilde{\sigma }}_{ik}={\tilde{\sigma }}_{i,k+1}$ or ${\tilde{\sigma }}_{jk}={\tilde{\sigma }}_{j,k+1}$. Thus, we mainly consider the cases of ${\tilde{\sigma }}_{jk}\ne {\tilde{\sigma }}_{j,k+1}$ and ${\tilde{\sigma }}_{ik}\ne {\tilde{\sigma }}_{i,k+1}$.

If ${\tilde{\sigma }}_{jk}>{\tilde{\sigma }}_{j,k-1}$, we have$$\begin{array}{c}\text{(B.5)}& {\tilde{\sigma }}_{jk}-{\tilde{\sigma }}_{j,k-\mathrm{1}}>\mathrm{0},\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{jk}^{\mathrm{2}}+{\widehat{S}}_u^{\mathrm{2}}}<\sqrt{\left({\tilde{\rho }}^{\mathrm{2}}-\mathrm{1}\right)r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j,k-\mathrm{1}}^{\mathrm{2}}+{\widehat{S}}_u^{\mathrm{2}}}.\end{array}$$It follows that ${\tilde{\sigma }}_{i,k+1}<{\tilde{\sigma }}_{ik}$ according to Eq. (B.3), and therefore ${\tilde{\sigma }}_{j,k+1}>{\tilde{\sigma }}_{jk}$ from Eq. (B.2). Similarly, if ${\tilde{\sigma }}_{jk}<{\tilde{\sigma }}_{j,k-1}$, we have ${\tilde{\sigma }}_{j,k+1}<{\tilde{\sigma }}_{jk}$. Hence, the sequence $\{{\tilde{\sigma }}_{jk}\}$ is monotonic. Similarly, it can be proved that the sequence $\{{\tilde{\sigma }}_{ik}\}$ is also monotonic.

According to Eqs. (22) and (24), we have ${\tilde{\sigma }}_{ik}\le \sqrt{{\widehat{S}}_u^{\mathrm{2}}}$ and ${\tilde{\sigma }}_{ik}\le \sqrt{{\widehat{S}}_v^{\mathrm{2}}}$. Note that, as standard deviation of a random variable, it is mandatory that ${\tilde{\sigma }}_{ik}\ge \mathrm{0}$ and ${\tilde{\sigma }}_{jk}\ge \mathrm{0}$. Thus, $\left\{{\tilde{\sigma }}_{ik},{\tilde{\sigma }}_{jk}\right\}$ are bounded. A bounded monotonic sequence should be convergent. Proposition 7 is proved.

Proof of Proposition 8.

Deduce the relationship between ${\tilde{\mu }}_{jk}$ and ${\tilde{\mu }}_{j,k-1}$, which yields$$\begin{array}{}\text{(B.6)}& {\tilde{\mu }}_{jk}={\stackrel{-}{x}}_v+\frac{\varphi \left(z_j\right)\sqrt{{\tilde{\sigma }}_{jk}^{\mathrm{2}}+r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{jk}}}{P_j}-r_{ij}{\tilde{\mu }}_{ik}={\stackrel{-}{x}}_v+\frac{\varphi \left(z_j\right)\sqrt{{\tilde{\sigma }}_{jk}^{\mathrm{2}}+r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{jk}}}{P_j}-r_{ij}\left({\stackrel{-}{x}}_u+\frac{\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_{ik}^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j,k-\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ji}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{jk}}}{P_i}-r_{ji}{\tilde{\mu }}_{j,k-\mathrm{1}}\right)={∆}_k+r_{ij}{r}_{ji}{\tilde{\mu }}_{j,k-\mathrm{1}},\end{array}$$where ${∆}_k=\left({\stackrel{-}{x}}_v-r_{ij}{\stackrel{-}{x}}_u\right)+\varphi \left(z_j\right)\sqrt{{\tilde{\sigma }}_{jk}^{\mathrm{2}}+r_{ij}^{\mathrm{2}}{\tilde{\sigma }}_{ik}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ij}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{jk}}/P_j-r_{ij}(\varphi \left(z_i\right)\sqrt{{\tilde{\sigma }}_{ik}^{\mathrm{2}}+r_{ji}^{\mathrm{2}}{\tilde{\sigma }}_{j,k-\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\tilde{\rho }r}_{ji}{\tilde{\sigma }}_{ik}{\tilde{\sigma }}_{jk}}/P_i)$. By iteration, it yields that ${\lim }_{k\to \infty }{\tilde{\mu }}_{jk}={\lim }_{k\to \infty }{∆}_k/(\mathrm{1}-r_{ij}{r}_{ji})$. As $\left\{{\tilde{\sigma }}_{ik},{\tilde{\sigma }}_{jk}\right\}$ are convergent, ${\lim }_{k\to \infty }{∆}_k$ exists so that $\left\{{\tilde{\mu }}_{ik}\right\}$ is convergent. Similarly, the sequence $\left\{{\tilde{\mu }}_{ik}\right\}$ is also convergent. This proves Proposition 8.

References

[1] R. Ernst, B. Kamrad, "Estimating demand by using sales information: inaccuracies encountered," European Journal of Operational Research, vol. 174 no. 2, pp. 675-688, DOI: 10.1016/j.ejor.2004.12.024, 2006.

[2] Y. Qin, R. Wang, A. J. Vakharia, Y. Chen, M. M. H. Seref, "The newsvendor problem: review and directions for future research," European Journal of Operational Research, vol. 213 no. 2, pp. 361-374, DOI: 10.1016/j.ejor.2010.11.024, 2011.

[3] J. Dhandapani, R. Uthayakumar, "Multi-item EOQ model for fresh fruits with preservation technology investment, time-varying holding cost, variable deterioration and shortages," Journal of Control and Decision, vol. 4 no. 2, pp. 70-80, 2017.

[4] M. Bijvank, I. F. A. Vis, "Lost-sales inventory theory: a review," European Journal of Operational Research, vol. 215 no. 1,DOI: 10.1016/j.ejor.2011.02.004, 2011.

[5] S. Nahmias, "Demand estimation in lost sales inventory systems," Naval Research Logistics (NRL), vol. 41 no. 6, pp. 739-757, DOI: 10.1002/1520-6750(199410)41:6<739::AID-NAV3220410605>3.0.CO;2-A, 1994.

[6] N. Agrawal, S. A. Smith, "Estimating negative binomial demand for retail inventory management with unobservable lost sales," Naval Research Logistics (NRL), vol. 43 no. 6, pp. 839-861, DOI: 10.1002/(SICI)1520-6750(199609)43:6<839::AID-NAV4>3.0.CO;2-5, 1996.

[7] S. Netessine, F. Zhang, "Positive vs. negative externalities in inventory management: Implications for supply chain design," Manufacturing and Service Operations Management, vol. 7 no. 1, pp. 58-73, DOI: 10.1287/msom.1040.0058, 2005.

[8] Z.-J. M. Shen, X. Su, "Customer behavior modeling in revenue management and auctions: a review and new research opportunities," Production Engineering Research and Development, vol. 16 no. 6, pp. 713-728, DOI: 10.3401/poms, 2007.

[9] S. S. Anand, J. G. Hughes, D. A. Bell, A. R. Patrick, "Tackling the cross-sales problem using data mining," Proceedings of PAKDD97, pp. 331-343, .

[10] R.-Q. Zhang, L.-K. Zhang, W.-H. Zhou, R. Saigal, H.-W. Wang, "The multi-item newsvendor model with cross-selling and the solution when demand is jointly normally distributed," European Journal of Operational Research, vol. 236 no. 1, pp. 147-159, DOI: 10.1016/j.ejor.2014.01.006, 2014.

[11] R. A. Fisher, "Properties and applications of Hh functions," Introduction to Mathematical Tables, vol. 1, pp. xxiv-xxxv, 1931.

[12] A. C. Cohen, "On Estimating the Mean and Standard Deviation of Truncated Normal Distributions," Journal of the American Statistical Association, vol. 44 no. 248, pp. 518-525, DOI: 10.1080/01621459.1949.10483324, 1949.

[13] A. K. Gupta, "Estimation of the mean and standard deviation of a normal population from a censored sample," Biometrika, vol. 39 no. 3-4, pp. 260-273, DOI: 10.1093/biomet/39.3-4.260, 1952.

[14] E. Berk, Ü. Gürler, R. A. Levine, "Bayesian demand updating in the lost sales newsvendor problem: a two-moment approximation," European Journal of Operational Research, vol. 182 no. 1, pp. 256-281, DOI: 10.1016/j.ejor.2006.08.035, 2007.

[15] H. S. Heese, J. M. Swaminathan, "Inventory and sales effort management under unobservable lost sales," European Journal of Operational Research, vol. 207 no. 3, pp. 1263-1268, DOI: 10.1016/j.ejor.2010.06.016, 2010.

[16] A. Jain, N. Rudi, T. Wang, "Demand estimation and ordering under censoring: Stock-out timing is (almost) all you need," Operations Research, vol. 63 no. 1, pp. 134-150, DOI: 10.1287/opre.2014.1326, 2015.

[17] A.-L. Beutel, S. Minner, "Safety stock planning under causal demand forecasting," International Journal of Production Economics, vol. 140 no. 2, pp. 637-645, DOI: 10.1016/j.ijpe.2011.04.017, 2012.

[18] A.-L. Sachs, S. Minner, "The data-driven newsvendor with censored demand observations," International Journal of Production Economics, vol. 149 no. 1, pp. 28-36, DOI: 10.1016/j.ijpe.2013.04.039, 2014.

[19] A. Azadeh, R. Arab, S. Behfard, "An adaptive intelligent algorithm for forecasting long term gasoline demand estimation: the cases of usa, canada, japan, kuwait and iran," Expert Systems with Applications, vol. 37 no. 12, pp. 7427-7437, DOI: 10.1016/j.eswa.2010.03.009, 2010.

[20] A. P. Dempster, N. M. Laird, D. B. Rubin, "Maximum likelihood from incomplete data via the em algorithm," Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 39 no. 1, 1977.

[21] R. Dahiya, R. Korwar, "Maximum likelihood estimates for a bivariate normal distribution with missing data," The Annals of Statistics, vol. 8 no. 3, pp. 687-692, 1980.

[22] K. Adamidis, S. Loukas, "Ml estimation in the bivariate poisson distribution in the presence of missing values via the em algorithm," Journal of Statistical Computation and Simulation, vol. 50 no. 3-4, pp. 163-172, DOI: 10.1080/00949659408811608, 1994.

[23] J. Huber, A. Gossmann, H. Stuckenschmidt, "Cluster-based hierarchical demand forecasting for perishable goods," Expert Systems with Applications, vol. 76 no. 15, pp. 140-151, DOI: 10.1016/j.eswa.2017.01.022, 2017.

[24] D. Besanko, S. Gupta, D. Jain, "Logit demand estimation under competitive pricing behavior: an equilibrium framework," Management Science, vol. 44 no. 11, pp. 1533-1547, DOI: 10.1287/mnsc.44.11.1533, 1998.

[25] R. Anupindi, M. Dada, S. Gupta, "Estimation of consumer demand with stock-out based substitution: An application to vending machine products," Marketing Science, vol. 17 no. 4, pp. 406-423, DOI: 10.1287/mksc.17.4.406, 1998.

[26] S. A. Smith, N. Agrawal, "Management of multi-item retail inventory systems with demand substitution," Operations Research, vol. 48 no. 1, pp. 50-64, DOI: 10.1287/opre.48.1.50.12443, 2000.

[27] A. G. Kök, M. L. Fisher, "Demand estimation and assortment optimization under substitution: methodology and application," Operations Research, vol. 55 no. 6, pp. 1001-1021, DOI: 10.1287/opre.1070.0409, 2007.

[28] C. T. Conlone, J. H. Mortimer, "Demand estimation under incomplete product availability," American Economic Journal: Microeconomics, vol. 5 no. 4, 2013.

[29] M. Wan, Y. Huang, L. Zhao, T. Deng, J. C. Fransoo, "Demand estimation under multi-store multi-product substitution in high density traditional retail," European Journal of Operational Research, vol. 266 no. 1, pp. 99-111, DOI: 10.1016/j.ejor.2017.09.014, 2018.

[30] R. Agrawal, T. Imielinski, A. Swami, "Mining association rules between sets of items in large databases," ACM SIGMOD Record, vol. 22 no. 2, pp. 207-216, 1993.

[31] R. Addo-Tenkorang, P. T. Helo, "Big data applications in operations/supply-chain management: a literature review," Computers & Industrial Engineering, vol. 101 no. 11, pp. 528-543, DOI: 10.1016/j.cie.2016.09.023, 2016.