When does skewness matter in robust inventory management?

Abstract

We investigate the impact of skewed demand on robust inventory management. Including skewness in calculations leads to cubic constraints that prevent a standard two-stage method from deriving an explicit and tractable objective function. To overcome this roadblock, we propose a joint optimization method to directly derive the robust solution in a closed form. The joint optimization method is widely applicable to various robust inventory models with different ambiguity sets. The notable advantage is that we can obtain the final solution without deriving the objective function, so many tedious intermediate steps are circumvented. We conduct numerical experiments on industry data to demonstrate that our moment-based policies deliver more consistent performance than divergence-based policies. After obtaining various closed-form solutions, we demonstrate that including skewness in the model improves the profit generated by the robust optimal order quantity if either demand is bounded or the cost-to-price ratio is low, even though the sample moments used may not equal the population moments. Furthermore, under these conditions in which skewness should be included in the calculations, the firm’s expected profit, under the most unfavorable distribution, increases with variance but decreases with skewness. This result stands in contrast to findings in the economics and finance literature, where distributional ambiguity is not considered.

1. Introduction

The vast majority of studies (and textbooks) of stochastic inventory management assume a known probability distribution for the demand faced by a firm. However, in many practical cases, the actual demand distribution is difficult to predict accurately (Lilien et al., 2017). Hence, some researchers have sought the ordering decision that maximizes the firm’s profit under the worst-case distribution (i.e., under the extreme or the most unfavorable distribution). This is known as the robust inventory management problem. Its earliest representative work is a seminal paper published by Scarf (1958). We study the robust inventory management problem for a profit-maximizing newsvendor.

Sales data from a large Chinese beverage company that offers 72 stock-keeping units (SKUs) provides an example of the robust inventory management problem because these data reveal uncertain demand distributions (https://www.coap.online/competitions/1). The data distributions for two representative SKUs (SKU17 and SKU23) over 181 days (from January 1st to June 30th, 2018) are shown in Figure 1. Readers can refer to the Excel file in the online E-Companion with the name SKUDataDisclosure.xlsx to find the details. The distribution patterns are completely different: the cumulative demand curve in Figure 1(a) initially accumulates rapidly and then increases gradually to $100\mathrm{\%}$ , whereas the curve in Figure 1(b) accumulates slowly at first (ignoring the mass on the left that represents days with zero demand) and then increases quickly to $100\mathrm{\%}$ . This makes accurately fitting them into any known distribution difficult. Consequently, the beverage company needs to solve the distribution-free (also called nonparametric) newsvendor inventory problem.

Figure 1.

Demands for SKU17 and SKU23 of the case company from January 1, 2018, to 30 June 30, 2018. (a) Sample mean $276.6$ , standard deviation $284.2$ , and skewness $1.616$ and (b) sample mean $\mathrm{70,982.8}$ , standard deviation $\mathrm{50,569.3}$ , and skewness $-0.378$ . SKU = stock-keeping unit.

An important feature of robust inventory models is their ambiguity set, which critically depends on what partial information is available (Long et al., 2024). When facing distributional ambiguity, firms usually make decisions based on a parsimonious set of input parameters, giving rise to a growing stream of literature on robust inventory management (see Lu and Shen, 2021, for an updated literature review). Among these parameters, the mean and variance, or equivalently the first two moments, are the two most commonly assumed known statistics (e.g., Fu et al., 2026, 2018; Gallego and Moon, 1993; Xing et al., 2022).

However, skewness, which is widely investigated in the finance literature (e.g., Chiu, 2005; Kumar et al., 2024; Lau et al., 1989) and measures the asymmetry of a distribution, thereby reflecting downside risk, is a factor that should be considered because the demand data from industries such as server computers, automotive, discrete manufacturing, and so on, usually exhibit left or right skew (e.g., Cheng et al., 2012; Natarajan et al., 2018; Saldanha et al., 2023). For example, for our beverage company, the demand for SKU17 is positively skewed and that of SKU23 is negatively skewed. The skewness makes the confidence interval (which is based on symmetry) no longer correct, resulting in asymmetric error probabilities. Specifically, if skewness is negative, the left tail of the demand is long and thin, but the right tail is short and fat; if the skewness is positive, the opposite occurs. Hence, skewness may generate unpredictable errors for inventory management, that is, ordering either more or less than optimal (Zhang et al., 2020). Nevertheless, when skewness is considered, a traditional two-stage method proposed by Scarf (1958) using known mean and variance encounters difficulty in finding a closed-form solution because of cubic constraints that make the analysis more challenging. This has caused the impact of skewness on robust inventory management to remain unexplored in the literature.

Motivated by the popularity of skewness in finance and economics and by the fundamental importance of considering distributionally robust inventory management in the newsvendor model, we address the following research questions: How can we explicitly formulate and efficiently solve a distribution-free newsvendor model if the first three moments are known? What would the most unfavorable distribution be, and how many units should be ordered to optimize the profits? When considering skewness, how do our proposed robust solutions perform in comparison to other models? We perform a systematic analysis based on a joint optimization approach to derive the optimal robust inventory policy after estimating the mean, variance, and skewness of observed data for which the distribution is unknown.

The rest of the paper is organized as follows. Subsection 1.1 reviews the related literature, and subsection 1.2 states our contributions. In section 2, we revisit the newsvendor model with unknown demand distribution, prove the feasibility of the joint optimization approach by applying duality theory, and present some useful conclusions. Section 3 specifies the worst-case distributions, the robust inventory decisions, and the corresponding optimized utility when considering mean, variance, and skewness under bounded demand. In section 4, we consider unbounded demand by investigating the effects of variance and skewness. Section 5 presents three extensions by applying the joint optimization process to two data-driven alternative models and to the conditional value at risk (CVaR) model. Section 6 examines the performance of the proposed robust policy using numerical experiments with both industry and simulated data. Section 7 concludes the paper.

1.1. Related literature

Our work belongs to the research stream that investigates the impact of higher moments on decisions. Markowitz (1952) made an initial attempt to develop this effect on portfolio management in economics and finance, inspiring numerous subsequent studies. The typical approach starts with a presumed demand distribution and then computes the higher moments of the profit or investment return. A consensus in this literature stream is that the mean pertains to the expected return, the variance pertains to volatility, and the skewness pertains to the risk—upside risk if positive, downside if negative (e.g., Briec et al., 2007; Samuelson, 1970; Theodossiou and Savva, 2016). Therefore, to be more profitable, the investor’s objective function allocates positive weights to the mean and skewness but a negative weight to the variance. Inspired by the economics and finance methodology, Zhang et al. (2020) modified the newsvendor’s objective function as a weighted average of the computed mean, variance, and skewness (see equation (39) in Zhang et al., 2020: 1411) by using exogenous weights (with a positive weight being allocated to skewness and a negative weight being allocated to variance). In contrast, we find that these weights should not have fixed signs. We identify conditions for which firms should endogenously allocate negative weight to skewness (so as to reduce the loss) but positive weight to variance.

The inventory literature also considers CVaR (e.g., Chen et al., 2009) or value at risk (e.g., Kouvelis et al., 2021) as alternative objectives. However, the presumed demand distribution critically affects the values of the computed moments. For instance, the computed skewness function has a V shape under a uniform distribution but an S shape under an exponential distribution (see Figures 1 and 2 in Zhang et al., 2020: 1404 and 1405). Obviously, these contrasting shapes of the skewness can produce disparate effects on inventory management. Moreover, if the demand distribution is unavailable or inaccurate, inventory managers may suffer significant losses due to severe model errors. Worse still, the differing shapes caused by skewness will further exacerbate these model errors. Therefore, distribution-free inventory management models have naturally received consistent praise from inventory scholars (e.g., Natarajan et al., 2018; Saldanha et al., 2023). As such, we assume that the first three moments of demand (rather than the entire distribution) have been calculated from collected data. In the analysis, we first endogenously identify the firm’s most unfavorable demand distribution using a joint optimization method and then derive the firm’s optimal robust inventory level. We also use this distribution as a reference to evaluate the performance of other policies.

Among the distribution-free research in inventory management, a significant achievement is the pioneering research of Scarf (1958), who formulates a two-stage model. The first stage computes the firm’s worst-case expected profit for any given inventory decision. Hereafter, we refer to the worst-case expected profit as the utility to highlight the use of the max–min decision rule and to keep the terminology concise. Given the expected profit function based on the worst-case distribution determined in the first stage, the second stage optimizes the inventory decision using the robust max–min decision rule. Because the analysis in the second stage involves the first-order condition (FOC), the firm’s utility function must be obtained explicitly in the first stage. To achieve this goal, the most popular method has been to reformulate the first-stage model as a deterministic model via duality (e.g., Fu et al., 2026; Li and Kirshner, 2021; Xing et al., 2022) or to apply the Cauchy inequality (Gallego and Moon, 1993).

However, closed-form solutions can be difficult to obtain if the information exceeds the mean and variance (Das et al., 2021) or if the firm’s utility becomes more complicated than that in the standard newsvendor model (e.g., Fu et al., 2018; Govindarajan et al., 2021; Han et al., 2014). First, in relation to our mean–variance–skewness robust model, the insurance and finance literature has studied the bounds on the expected linear loss (equivalently, the expected lost sales quantity) based on the first three moments. To be specific, Jansen et al. (1986) and Schepper and Heijnen (2007) provided the most accurate and relevant expressions. Nonetheless, as we restate and demonstrate these known bounds in Appendix A of the online E-Companion, the firm’s objective function is overly complex in some cases, making the FOC intractable. Second, the issue of intractable FOCs also exists in a similar model in which the mean and the $t$ -th moment are known. We refer to this model as the $1+t$ model. If we follow the traditional two-stage method, the FOC is intractable. For example, Das et al. (2021: 1081) summarized the challenge pertaining to the $1+t$ model as the following: “In general, for these problems, there is an absence of closed-form solutions and hence finding an explicit representation of the distribution $F^{*}$ does not appear to be straightforward.”

To overcome the problems caused by higher moments and the challenge of deriving the FOCs, we propose a joint optimization approach to successfully characterize the robust inventory solutions of a newsvendor and to facilitate the subsequent analysis with known mean, variance, and skewness. This joint optimization method can be widely applied to many robust inventory models as long as we appropriately change the ambiguity set. For example, we can solve the formidable $1+t$ model because we directly attack the final answer instead of going through tedious intermediate steps (which merely produce an implicit objective function). Our case study also shows that our robust solution performs well even when the input data of moments is based on a partial (rather than the entire) history.

Finally, we are unaware of any published article that unambiguously and distinctly applies joint optimization to derive analytical solutions for moment-based problems. To our best knowledge, the literature still follows a two-stage approach. Additionally, many articles jointly optimize complementary or substitutable objectives in which they intend to design coordinated actions for the firm. For example, Bensoussan et al. (2019) consider price–inventory and Sloan and Shanthikumar (2000) consider production–maintenance joint optimization challenges. On the contrary, the approach we propose directly discards the above two-stage method and transforms it into a single-stage problem that is easier to solve through joint optimization; specifically, our method determines the worst-case distribution and the corresponding order quantity within a single stage. This enables us to clearly obtain the extreme distributions as well as a vector of actions (i.e., the optimal order quantity). On the other hand, from the perspective of game theory, the above two-stage problem can also be regarded as a zero-sum game in which the firm chooses an order quantity to maximize profits, while the opponent of the firm (i.e., nature) chooses a distribution to reduce the firm’s profit. Our proposed joint optimization approach simplifies the problem by transforming the original intricate two-player game between nature and the firm into a single-player’s (the firm’s) decision problem, thereby making the model easier to solve.

1.2. Our contributions

We highlight our contributions and results as follows.

We develop a robust newsvendor model based on the first three moments. We demonstrate that skewness is useful for describing asymmetric demand distributions if demand is bounded or if the cost-to-price ratio is sufficiently high. Otherwise—unbounded demand with a low or moderate cost-to-price ratio—skewness does not contribute to specifying the demand distribution, and the robust production quantity should follow Scarf’s (1958) rule. This managerial insight should assist firms in deciding whether to include skewness in their robust inventory management process.

We show that if demand is bounded or if the cost-to-price ratio is sufficiently low, then firms should seek larger variance and smaller skewness when calculating optimal order quantities. This result contradicts standard results in economics and finance but arises when the external demand follows the most unfavorable distribution or a beta distribution.

We propose a joint optimization method to derive closed-form solutions. When the third moment is included in the robust inventory model, we show that the traditional method encounters a roadblock in solving the first-stage problem (in which the firm needs to determine the worst-case objective function): of the five cases produced by this method, two require finding the root of an endogenous cubic equation. As a result, the FOC of the firm’s utility function is unavailable in explicit form, hindering the subsequent analysis. By contrast, we propose a joint optimization method that bypasses the first-stage problem to efficiently and effectively derive order quantities for various cases in closed form. We also identify the most unfavorable distribution, which serves as a useful benchmark of performance.

We use industrial data to demonstrate the effectiveness of our robust policies. We conduct an applied numerical experiment using a data set from a large beverage company in China. We demonstrate that after including skewness in the planning process, the moment-based policy outperforms rival policies, including a policy based on an empirical distribution and the Kullback–Leibler (KL) divergence.

2. Model and methodology

We consider a single-period newsvendor model in which the external demand that the firm receives is represented by a nonnegative random variable $\tilde{\theta }$ . Hereafter, we use a tilde to indicate a random variable and notation without a tilde to indicate a realized value. The firm’s selling price is $p$ and the production cost is $c$ , in which $c<p$ . The firm produces $q$ units, and then the external demand arrives. The realized demand above the available inventory is lost, and any unsold inventory is salvaged to zero. We define the ex-post profit as $Z(\theta ,q)=pmin(\theta ,q)-cq$ . Let $F$ be the unknown cumulative distribution function of $\tilde{\theta }$ .

Although unbounded demand is common in moment-based inventory models, it can be inconsistent with many practical circumstances, resulting in overstocking problems identified in the literature. Thus, we first consider bounded demand (section 3) and then address the issue of unbounded demand in section 4. We scale the bounded demand to the closed unit interval $[0,1]$ to streamline the notation. Although the exact form of $F$ remains unknown, we assume that the first four moments ( ${\int }_0^1{\theta }^{i}dF(\theta )=m_i$ , where $i\in \{0,1,2,3\}$ ) are known. By default, $m_0=1$ represents the total probability, $m_1=\mu$ represents the mean, $m_2={\mu }^2+{\sigma }^2$ , where ${\sigma }^2$ represents the variance, and $m_3$ determines the skewness. The values of the $m_i$ represent the partial information available to the firm. Let ${\mathfrak{R}}_{+}$ be the collection of nonnegative real numbers and $\mathcal{M}({\mathfrak{R}}_{+})$ be the nonnegative regular Borel measures. Based on the known moments, we define the following ambiguity set:

$$\begin{array}{r}\mathrm{\Omega }=\{F(\cdot )\in \mathcal{M}\left({\mathfrak{R}}_{+}\right)|{\int }_0^1{\theta }^{i}dF\left(\theta \right)=m_i,i=0,1,2,3\}\end{array}$$

(1)

Throughout this paper, we let $\mathbb{E}$ be the expectation operator. We define skewness $S_{\mathrm{k}}$ as follows:

$$\begin{array}{rl}{S}_{\mathrm{k}}& =\mathbb{E}\left[{\left(\frac{\tilde{\theta }-\mu }{\sigma }\right)}^3\right]=\frac{m_3-3\mu {\sigma }^2-{\mu }^3}{{\sigma }^3}\\ & =\frac{m_3-3m_1{m}_2+2m_1^3}{{(m_2-m_1^2)}^{3/2}}\end{array}$$

(2)

If $S_{\mathrm{k}}$ is positive (resp., negative), then the demand distribution is positively (resp., negatively) skewed.

A critical issue is to determine whether the set $\mathrm{\Omega }$ in equation (1) contains an infinite number of feasible distributions. Lemma 1 answers this question by providing bounds on each moment based on its lower moments.

Lemma 1 Jansen et al., 1986

When $\tilde{\theta }$ has its range within the unit interval, the moments must satisfy (i) $0\le m_1\le 1$ ; (ii) $m_1^2\le m_2\le m_1$ ; and (iii) $m_2^2/m_1\le m_3\le (m_2-m_2^2+m_1{m}_2-m_1^2)/(1-m_1)$ .

Whenever one of the inequalities in Lemma 1 holds with the equal sign, the ambiguity set $\mathrm{\Omega }$ in equation (1) contains only one feasible distribution. Readers may refer to Remark EC.1 in the online E-Companion for details. In such cases with only one feasible distribution, the analysis is reduced to the unambiguous model. Thus, to avoid uninteresting cases, we hereafter assume that Lemma 1 strictly holds. The firm’s objective is to solve the distributional robust optimization model

$$Z^{*}=\underset{q}{max}\underset{F\in \mathrm{\Omega }}{inf}\left\{{\int }_0^{1}Z(\theta ,q)dF\left(\theta \right)\right\}$$

(3)

by finding $q$ and $F$ .

2.1. Duality and joint optimization

For clarity, we define a few identities. Let

$$D=\underset{F\in \mathrm{\Omega }}{inf}{\int }_0^{1}Z(\theta ,q)dF\left(\theta \right)\stackrel{\mathrm{def}}{=}{Z}_{wst}\left(q\right)$$

be the dual in measure space and $Z_{wst}(q)$ be the firm’s utility or worst-case expected profit. Although primal and dual can be used interchangeably, the literature (e.g., Hettich and Kortanek, 1993) established a convention that the model with an infinite number of decision variables and a finite number of constraints is called the dual; whereas the model with an infinite number of constraints and a finite number of decision variables is called the primal. The traditional two-stage method for solving equation (3) first derives $Z_{wst}(q)$ and then solves the FOC $\mathrm{\partial }{Z}_{wst}(q)/\mathrm{\partial }q=0$ . Let $F_{wst}(\cdot |q)$ be the extreme distribution that attains the infimum for a given $q$ , that is, $Z_{wst}(q)={\int }_0^{1}Z(\theta ,q)d{F}_{wst}(\theta |q)$ . Let $q^{*}$ be the order quantity maximizing $Z_{wst}(q)$ . We call $q^{*}$ the robust optimal solution. Correspondingly, we refer to $F^{*}(\cdot )=F_{wst}(\cdot |q^{*})$ as the firm’s most unfavorable distribution, so $Z^{*}=Z_{wst}(q^{*})={\int }_0^{1}Z(\theta ,q^{*})d{F}^{*}(\theta )$ is the firm’s optimized utility.

The literature also considers the dual in finite sequence space. The sequence $\{\lambda (\theta )\}$ is said to be a general finite sequence if $\lambda (\theta )\ge 0$ for any $\theta \in [0,1]$ , but only a finite number of elements of this sequence can be positive (Hettich and Kortanek, 1993: 386, 389, 390). Hence, we can apply summation to define the following ambiguity set:

$${\mathrm{\Omega }}^{\mathrm{\prime }}=\left\{\lambda \left(\theta \right)|\sum _{\theta }{\theta }^i\lambda \left(\theta \right)=m_i,i=0,1,2,3;\theta \in [0,1]\right\}$$

It must hold that ${\mathrm{\Omega }}^{\mathrm{\prime }}\subset \mathrm{\Omega }$ since ${\mathrm{\Omega }}^{\mathrm{\prime }}$ includes only discrete distributions while the set $\mathrm{\Omega }$ includes discrete, mixed, and continuous distributions. The dual in general finite sequence space equals

$$D^{\mathrm{\prime }}=\underset{\lambda \left(\theta \right)\in {\mathrm{\Omega }}^{\mathrm{\prime }}}{min}\sum _{\theta }Z(\theta ,q)\lambda \left(\theta \right)$$

where $D\le D^{\mathrm{\prime }}$ must hold. A sufficient condition to establish $D=D^{\mathrm{\prime }}$ is that the ex-post payoff function $Z(\theta ,q)$ is continuous. In the newsvendor model, $Z(\theta ,q)$ $=pmin(\theta ,q)-cq$ is continuous in $\theta$ , resulting in $D=D^{\mathrm{\prime }}$ .

For any given $q$ , we formulate the following semi-infinite programming (SIP) model $P$ :

$$\begin{array}{rl}P=& \underset{\overrightarrow{y}}{max}\{y_0+y_1{m}_1+y_2{m}_2+y_3{m}_3\}\\ \text{s.t. }& y_0+y_1\theta +y_2{\theta }^2+y_3{\theta }^3\le Z(\theta ,q),\mathrm{\forall }\theta \in [0,1],\end{array}$$

where the decision variables are $y_i$ , $i=0,1,2,3$ . As each $\theta$ defines one constraint and $\theta$ can be any real number within the unit interval, the number of constraints is infinite. The technical convenience is that model $D^{\mathrm{\prime }}$ (rather than model $D$ ) and model $P$ form a dual-primal pair. The primal model matching with model $D$ is inconvenient for subsequent analysis. Weak duality $P\le D^{\mathrm{\prime }}$ holds, but we now establish $D=D^{\mathrm{\prime }}=P$ .

Lemma 2

If (i) $Z(\theta ,q)$ is continuous with respect to $\theta$ for any given $q$ , (ii) $Z(\theta ,q)$ is finite for any given $(\theta ,q)$ , and (iii) the sequence $\{1,m_1,m_2,m_3\}$ strictly satisfies Lemma 1, then $D=D^{\mathrm{\prime }}=P$ .

The three joint sufficient conditions in Lemma 2 are mild. The first condition ensures that $D=D^{\mathrm{\prime }}$ (i.e., the infima in measure space and in finite sequence space are equal). The second and third conditions ensure that $D^{\mathrm{\prime }}=P$ (i.e., duality in the finite sequence). These conditions hold in many supply chain problems. For example, due to lost sales or limited liability, the ex-post order quantity is finite, making the ex-post payoff finite. ¹

When $D=D^{\mathrm{\prime }}=P$ holds, we can reformulate the firm’s model (3) as follows:

$$\begin{array}{rl}Z& =\underset{q}{max}\underset{y_i}{max}\left\{\sum _{i=0}^3{y}_i{m}_i\right\}=\underset{q,y_i}{max}\left\{\sum _{i=0}^3{y}_i{m}_i\right\}\end{array}$$

(4)

$$\begin{array}{rl}& \text{s.t. }{y}_0+y_1\theta +y_2{\theta }^2+y_3{\theta }^3\le Z(\theta ,q),\mathrm{\forall }\theta \in [0,1]\end{array}$$

(5)

The objective function in equation (4) is linear. Because the newsvendor model yields an ex-post payoff $Z(\theta ,q)$ that is continuous, concave, and finite with respect to $\theta$ , the feasible region in equation (5) is compact and convex, making the Karush–Kuhn–Tucker (KKT) conditions sufficient and necessary. By including $q$ in the KKT conditions, we surprisingly find that many roadblocks stemming from the traditional two-stage methods vanish. ²

The joint optimization method closely relates to zero-sum games. Any max–min model inherently represents a zero-sum game between the firm and nature. While the firm chooses actions to maximize the expected profit, nature chooses a distribution to minimize the expected profit. The available partial information defines nature’s action space and, hence, more available information reduces nature’s action space and improves the firm’s equilibrium utility. The extreme distribution $F_{wst}(\cdot |q)$ is nature’s best response to any given $q$ , whereas $F^{*}(\cdot |q^{*})$ is nature’s equilibrium strategy (or the most unfavorable demand distribution). A zero-sum game perspective facilitates the analysis because in any equilibrium only a dominating strategy will be played. As the firm’s FOC characterizes the firm’s dominating strategy, the joint optimization method enables us to solve the equilibrium from a single player’s perspective. This is the main reason why our joint optimization method can bypass many intermediate steps and can derive analytical solutions that are otherwise unattainable by the traditional two-stage method.

3. Analytical results: Bounded demand

This section finds the optimal solutions and the most unfavorable distributions for four cases that span the problem’s domain. To facilitate the subsequent analysis, we define four useful thresholds $\{r_0,r_2^{\mathrm{\prime }},r_2,r_1\}$ :

$$\begin{array}{rl}{r}_0& =\frac{m_1{m}_3-m_2^2}{m_1+m_3-2m_2},r_2^{\mathrm{\prime }}=\frac{m_1^2}{m_2},\\ r_2& =\frac{m_2^2-m_1(m_2+m_3-m_1)}{m_2-m_3}\end{array}$$

(6)

$$\begin{array}{rl}{r}_1& =\{\begin{array}{ll}\frac{1}{2}+\frac{1}{2}\sqrt{1-\frac{4}{4+S_{\mathrm{k}}^2}},& \text{if }{S}_{\mathrm{k}}\le 0\\ \frac{1}{2}-\frac{1}{2}\sqrt{1-\frac{4}{4+S_{\mathrm{k}}^2}},& \text{if }{S}_{\mathrm{k}}\ge 0.\end{array}\end{array}$$

(7)

We note that $r_2^{\mathrm{\prime }}$ is the threshold in Scarf’s model to recommend zero quantity, whereas the other three thresholds involve information about skewness. This section’s results demonstrate why these four thresholds are crucial in our analysis.

Lemma 3

If the conditions of Lemma 1 strictly hold, then the four thresholds satisfy $0\le r_0\le r_1\le r_2^{\mathrm{\prime }}\le r_2\le 1$ .

The joint optimization model in equations (4) and (5) yields four cases for bounded demands, depending on the cost-to-price ratio’s relation to $r_0$ , $r_1$ , and $r_2$ . Each case is described in one of the following four propositions, which we present in ascending order of the cost-to-price ratio $c/p$ . First, Proposition 1 shows that if the cost-to-price ratio $c/p$ is sufficiently small, then the penalty for overstocking is so low that even the most unfavorable distribution cannot deter the firm from ordering the maximum quantity. The optimal order quantity decreases as the cost-to-price ratio increases in the subsequent propositions until, in Proposition 4, the cost-to-price ratio is so high that $q^{*}=0$ .

Proposition 1

If $c/p\le r_0$ , then the optimal solution satisfies $q^{*}=1$ , the optimized utility equals $Z^{*}=p{m}_1-c$ , and nature’s equilibrium strategy $F^{*}$ satisfies

$$F^{*}={\hat{F}}_1=\{\begin{array}{l}Pr(\tilde{\theta }=0)=1-r_2\\ Pr(\tilde{\theta }=\frac{m_2-m_3}{m_1-m_2})=r_2-r_0\\ Pr(\tilde{\theta }=1)=r_0.\end{array}$$

(8)

The second case involves the cubic equation

$$A(x)=\frac{c}{p}{m}_1{x}^3-\frac{2c}{p}{m}_2{x}^2+\frac{c}{p}{m}_{3}x+m_2^2-m_1{m}_3=0,$$

(9)

which depends on the cost-to-price ratio. Let ${\theta }_3$ be the largest real root of equation (9). Next, we use the known value of ${\theta }_3$ to specify another constant: ${\theta }_2=(m_3-{\theta }_3{m}_2)/(m_2-{\theta }_3{m}_1)$ .

Proposition 2

If $r_0<c/p\le r_1$ , then nature’s equilibrium strategy $F^{*}$ satisfies

$$F^{*}={\hat{F}}_2=\{\begin{array}{l}Pr(\tilde{\theta }={\theta }_1=0)=\frac{p-c}{p}-\frac{m_1-(c/p){\theta }_3}{{\theta }_2}\\ Pr(\tilde{\theta }={\theta }_2)=\frac{m_1-(c/p){\theta }_3}{{\theta }_2}\\ Pr(\tilde{\theta }={\theta }_3)=\frac{c}{p}.\end{array}$$

(10)

The firm’s optimal solution $q^{*}$ satisfies

$$q^{*}=\frac{2{\theta }_3^2}{3{\theta }_3-{\theta }_2}$$

(11)

and the optimized utility equals $Z^{*}=p{m}_1-c{\theta }_3$ .

In the first and second cases (Propositions 1 and 2), nature’s equilibrium strategies are three-point distributions with zero being the common realized value. While nature’s strategy in the second case depends on the cost-to-price ratio, that of the first case does not.

The third case involves a different cubic equation:

$$B(x)=B_3{x}^3+B_2{x}^2+B_{1}x+B_0=0,$$

where the four coefficients satisfy

$$\{\begin{array}{l}{B}_3=\frac{p-c}{p}(m_1-1)\\ B_2=\frac{p-c}{p}(1+m_1-2m_2)\\ B_1=\frac{p-c}{p}(m_2+m_3-2m_1)\\ B_0=\frac{c}{p}(m_3-m_2)-m_1{m}_2-m_1{m}_3+m_1^2+m_2^2.\end{array}$$

We must define roots ${\theta }_i$ , $i=4,5$ , and probability masses ${\lambda }_i$ , $i=4,5,6$ , before stating Proposition 3. Let ${\theta }_4$ be the smallest positive real solution of the cubic equation $B(x)=0$ . We define ${\lambda }_4=(p-c)/p$ as a relevant probability mass. We use both ${\lambda }_4$ and ${\theta }_4$ as inputs to specify another probability mass ${\lambda }_6$ :

$${\lambda }_6=\frac{{(m_1-{\lambda }_4{\theta }_4)}^2-(1-{\lambda }_4)(m_2-{\lambda }_4{\theta }_4^2)}{2m_1-2{\lambda }_4{\theta }_4+{\lambda }_4{\theta }_4^2-m_2+{\lambda }_4-1}.$$

Next, we use the known values of $({\lambda }_6,{\lambda }_4,{\theta }_4)$ as inputs to specify the two remaining constants:

$${\lambda }_5=\frac{c}{p}-{\lambda }_6\text{and}{\theta }_5=\frac{m_2-{\lambda }_4{\theta }_4^2-{\lambda }_6}{m_1-{\lambda }_4{\theta }_4-{\lambda }_6}.$$

We now present the third case.

Proposition 3

If $r_1<c/p\le r_2$ , then nature’s equilibrium strategy is the following three-point distribution:

$$F^{*}={\hat{F}}_3=\{\begin{array}{l}Pr(\tilde{\theta }={\theta }_4)=\frac{p-c}{p}={\lambda }_4\\ Pr(\tilde{\theta }={\theta }_5)=\frac{c}{p}-{\lambda }_6\\ Pr(\tilde{\theta }=1)={\lambda }_6.\end{array}$$

(12)

The firm’s robust optimal solution satisfies

$$q^{*}=\frac{{\theta }_4+{\theta }_5-2{\theta }_4^2}{2-3{\theta }_4+{\theta }_5},$$

(13)

and its optimized utility equals $Z^{*}=(p-c){\theta }_4$ .

Two observations can be made regarding Proposition 3. First, nature’s equilibrium strategy $F^{*}$ displays two characteristics: (i) a realized value $\theta =1$ (i.e., the upper bound of the demand), and (ii) the firm has a stockout with probability $c/p$ . Second, the information about skewness makes the firm less likely to cease production than in the mean–variance model. To be specific, recall that when $r_2^{\mathrm{\prime }}\le c/p$ , Scarf’s (1958) rule recommends a zero order quantity. In contrast, Proposition 3 shows that $q^{*}$ is still positive for $r_2\ge c/p\ge r_2^{\mathrm{\prime }}$ (since the thresholds satisfy $r_2\ge r_2^{\mathrm{\prime }}\ge r_1$ ).

Proposition 4

If $c/p>r_2$ , then the firm’s optimal solution satisfies $q^{*}=0$ , nature’s equilibrium strategy satisfies $F^{*}={\hat{F}}_1$ , and the optimized utility is $Z^{*}=0$ .

Nature plays the same equilibrium strategy ${\hat{F}}_1$ for both Propositions 1 and 4, but the firm chooses different order quantities. The optimal $q^{*}$ and $Z^{*}$ for each case with bounded demand using our scheme are summarized in Table 1.

Table 1.

Summary of results for Propositions 1 through 4.

$\mathbf{c}\mathbf{/}\mathbf{p}$	$\mathbf{(}\mathbf{0}\mathbf{,}{\mathbf{r}}_{\mathbf{0}}\mathbf{]}$	$\mathbf{(}{\mathbf{r}}_{\mathbf{0}}\mathbf{,}{\mathbf{r}}_{\mathbf{1}}\mathbf{]}$	$\mathbf{(}{\mathbf{r}}_{\mathbf{1}}\mathbf{,}{\mathbf{r}}_{\mathbf{2}}\mathbf{]}$	$\mathbf{[}{\mathbf{r}}_{\mathbf{2}}\mathbf{,}\mathbf{1}\mathbf{]}$
${\mathbf{q}}^{*}$	$1$	$\frac{2{\theta }_3^2}{3{\theta }_3-{\theta }_2}$	$\frac{{\theta }_4+{\theta }_5-2{\theta }_4^2}{2-3{\theta }_4+{\theta }_5}$	$0$
${\mathbf{Z}}^{*}$	$p{m}_1-c$	$p{m}_1-c{\theta }_3$	$(p-c){\theta }_4$	$0$

Example 1

If $m_1=1/2$ , $m_2=1/3$ , and $m_3=1/4$ (which are based on the standard uniform distribution), then the three relevant thresholds include $r_0=1/6$ , $r_1=1/2$ , and $r_2=5/6$ . Hence,

1. (i)

If $c/p\le 1/6$ , then $q^{*}=1$ ;

(ii)

If $1/6<c/p\le 1/2$ , then the firm’s robust production quantity follows equation (11);

(iii)

If $1/2<c/p\le 5/6$ , then the firm’s robust production quantity follows equation (13);

(iv)

If $c/p>5/6$ , then $q^{*}=0$ .

Example 1 illustrates the four cases shown in Propositions 1 to 4. In contrast, under Scarf’s rule, the firm ceases operations ( $q^{*}=0$ ) if $c/p\ge m_1^2/m_2=3/4$ .

4. Comparison analysis: Unbounded demand

We now assume that the firm’s demand is unbounded so that we can compare our model to Scarf’s (1958) model and to the semivariance model studied by Natarajan et al. (2018). We do this by modifying the ambiguity set as follows:

$$\begin{array}{rl}{\mathrm{\Omega }}_{\mathrm{\infty }}& =\{F(\cdot )\in \mathcal{M}\left({\mathfrak{R}}_{+}\right)|{\int }_0^{\mathrm{\infty }}{\theta }^{i}dF\left(\theta \right)=m_i,\\ & \phantom{{\int }_0^{\mathrm{\infty }}}i=0,1,2,3\},\end{array}$$

(14)

where subscript $\mathrm{\infty }$ indicates that demand is unbounded. The Stieltjes moment problem seeks sufficient and necessary conditions for $m_i$ such that the set ${\mathrm{\Omega }}_{\mathrm{\infty }}$ in equation (14) is nonempty. For the first four moments, the Stieltjes requirements include

$$\begin{array}{rl}& det\left(\begin{array}{cc}{m}_0& m_1\\ m_1& m_2\end{array}\right)=m_0{m}_2-m_1^2>0\text{and }\\ & det\left(\begin{array}{cc}{m}_1& m_2\\ m_2& m_3\end{array}\right)=m_1{m}_3-m_2^2>0.\end{array}$$

(15)

Hence, both $m_0{m}_2-m_1^2={\mu }^2+{\sigma }^2-{\mu }^2={\sigma }^2>0$ and $m_1{m}_3-m_2^2=\mu m_3-({\mu }^2+{\sigma }^2{)}^2>0$ must hold.

A noteworthy technical issue is that the set ${\mathrm{\Omega }}_{\mathrm{\infty }}$ in equation (14) is a weakly closed set. A common tool to make the ambiguity set closed in topology is to make the last condition an inequality. Readers can refer to Theorems 25.11 and 25.12 in Billingsley (1995) for more discussion about making the set closed. From the zero-sum game perspective, nature’s action space needs to be closed. As such, we define a larger set

$${\mathrm{\Omega }}_{\mathrm{\infty }}^{\mathrm{\prime }}=\{F(\cdot )\in \mathcal{M}\left({\mathfrak{R}}_{+}\right)|\begin{array}{l}{\int }_0^{\mathrm{\infty }}{\theta }^{i}dF\left(\theta \right)=m_i,i=0,1,2\\ \\ {\int }_0^{\mathrm{\infty }}{\theta }^{3}dF\left(\theta \right)\le m_3\end{array}\},$$

which satisfies ${\mathrm{\Omega }}_{\mathrm{\infty }}^{\mathrm{\prime }}\supset {\mathrm{\Omega }}_{\mathrm{\infty }}$ .

4.1. Skewness versus variance

We write the unbounded version of the joint optimization model as follows:

$$\begin{array}{rl}Z=& \underset{q,y_i}{max}\{y_0+y_1{m}_1+y_2{m}_2+y_3{m}_3\}\\ \text{s.t. }& y_0+y_1\theta +y_2{\theta }^2+y_3{\theta }^3\le Z(\theta ,q),\mathrm{\forall }\theta \ge 0.\end{array}$$

(16)

Naturally, $y_3\le 0$ must hold regardless of the sign of $y_1$ or $y_2$ (thus, only $y_3\le 0$ is required); otherwise, the SIP constraints of equation (16) fail when $\theta$ approaches positive infinity. Hence, the new constraint $y_3\le 0$ is automatically embedded as part of the joint optimization model in equation (16). From a primal-dual perspective, the constraint $y_3\le 0$ implies that ${\int }_0^{\mathrm{\infty }}{\theta }^{3}dF(\theta )\le m_3$ , where we change the $=$ sign to $\le$ . Thus, with unbounded demands, we can expand the ambiguity set from ${\mathrm{\Omega }}_{\mathrm{\infty }}$ to ${\mathrm{\Omega }}_{\mathrm{\infty }}^{\mathrm{\prime }}$ without altering the optimal solution. In contrast, with bounded demands, set $\mathrm{\Omega }$ in equation (1 ) is known to be closed in topology, and this technical issue vanishes. Hence, we use equations for all the moment conditions when demand is bounded.

Using set ${\mathrm{\Omega }}_{\mathrm{\infty }}^{\mathrm{\prime }}$ offers a few advantages. First, because set ${\mathrm{\Omega }}_{\mathrm{\infty }}^{\mathrm{\prime }}$ is closed in topology, it includes the two-point distributions ${\hat{F}}_4$ , given in equation (17) and ${\hat{F}}_5$ in equation (20), rather than their “converging and limiting versions.” The infimum becomes attained rather than asymptotically approached. Second, we can include $y_3\le 0$ as the new constraint for equation (16) and capitalize on the corresponding Lagrangian multiplier (see the proof of Proposition 5).

The robust optimal solution for the unbounded version requires two cases, as stated in the following propositions.

Proposition 5 When Skewness is Unimportant

(i) If $r_1\le c/p\le r_2^{\mathrm{\prime }}$ , then the firm’s most unfavorable distribution $F^{*}$ satisfies

$$F^{*}={\hat{F}}_4=\{\begin{array}{l}Pr(\tilde{\theta }=m_1-\sqrt{\frac{c}{p-c}}\sqrt{m_2-m_1^2})=\frac{p-c}{p}\\ Pr(\tilde{\theta }=m_1+\sqrt{\frac{p-c}{c}}\sqrt{m_2-m_1^2})=\frac{c}{p}.\end{array}$$

(17)

The firm’s robust optimal solution satisfies

$$q_{\mathrm{\infty }}^{*}=m_1+\frac{1}{2}\sqrt{m_2-m_1^2}(\sqrt{\frac{p-c}{c}}-\sqrt{\frac{c}{p-c}}),$$

(18)

and the optimized utility equals

$$Z_{\mathrm{\infty }}^{*}=(p-c)m_1-\sqrt{c(p-c)(m_2-m_1^2)}.$$

(19)

(ii) If $c/p>r_2^{\mathrm{\prime }}$ , then the firm’s most unfavorable distribution $F^{*}$ satisfies

$$F^{*}={\hat{F}}_5=\{\begin{array}{l}Pr(\tilde{\theta }=0)=\frac{m_2-m_1^2}{m_2}\\ Pr(\tilde{\theta }=\frac{m_2}{m_1})=\frac{m_1^2}{m_2}.\end{array}$$

(20)

The firm’s robust optimal solution satisfies $q_{\mathrm{\infty }}^{*}=0$ , and the optimized utility equals $Z_{\mathrm{\infty }}^{*}=0$ .

Proposition 5 shows that if $r_1\le c/p$ , then the firm should ignore the information about skewness. We now explain the implications of the threshold $r_1$ and the function $G(x)=(1-2x)/(\sqrt{x(1-x)})$ , which pertains to the skewness of the Bernoulli distribution (see the proof of Lemma 3). By using standardization (i.e., let $\tilde{\xi }=(\tilde{\theta }-m_1)/(\sqrt{m_2-m_1^2})$ such that $E(\tilde{\xi })=0$ and $\mathrm{Var}(\tilde{\xi })=1$ ), the two-point distribution ${\hat{F}}_4$ in equation (17) becomes

$${\hat{F}}_4=\{\begin{array}{l}Pr(\tilde{\xi }=-\sqrt{\frac{1-x}{x}})=1-x\\ Pr(\tilde{\xi }=+\sqrt{\frac{x}{1-x}})=x,\end{array}$$

where $x=c/p$ is the cost-to-price ratio. We notice that $\tilde{\xi }$ is a transformed Bernoulli distribution (with endogenous realizations determined by $x$ ), for which the skewness equals

$$G\left(x\right)=\frac{1-2x}{\sqrt{x(1-x)}}.$$

When the cost-to-price ratio is within the range that Proposition 5 specifies, the skewness constraint is nonbinding but the embedded constraint $y_3=0$ is binding. Nature’s equilibrium strategy ${\hat{F}}_4$ or ${\hat{F}}_5$ is identical to the corresponding case in Scarf’s model.

Proposition 6 When Skewness Matters

If $c/p<r_1$ , then the firm’s most unfavorable distribution is the three-point distribution characterized by Proposition 2. Both the firm’s robust optimal solution $q_{\mathrm{\infty }}^{*}$ and optimized utility $Z_{\mathrm{\infty }}^{*}$ follow Proposition 2.

When $c/p<r_1$ , the skewness constraint is binding, making the three-point distribution ${\hat{F}}_2$ in equation (10) emerge as nature’s equilibrium strategy. ³ A binding skewness constraint implies that the firm benefits from acquiring information about skewness. By using the implicit function theorem, we obtain the following comparative statics.

Corollary 1

If $c/p<r_1$ , then $\mathrm{\partial }{Z}^{*}/\mathrm{\partial }{m}_2>0$ and $\mathrm{\partial }{Z}^{*}/\mathrm{\partial }{m}_3<0$ .

Corollary 1 identifies that when $c/p<r_1$ , operations managers’ utility is increasing in variance but decreasing in skewness. This result stands in contrast to the financial literature (e.g., Briec et al., 2007) in which a positively skewed distribution of returns implies that an investor may face frequent small losses but occasional large gains. Since frequent small losses are manageable, while the occasional large gains are highly desirable, financial managers may assign positive weight to skewness to reflect their preference to limit the downside risk and to embrace upside potential. Alternatively, a negatively skewed demand distribution suggests that an operations manager may encounter frequent small stockouts but occasional large overstocks. If overstocking is significantly less costly than understocking (i.e., if the cost-to-price ratio is low), then operations managers can optimize the trade-off between expected understock and overstock costs by increasing the inventory level. Thus, operations managers may assign negative weight to skewness to indicate their preference when considering the worst-case demand distribution. In Appendix C of the online E-Companion, we provide additional numerical examples to validate Corollary 1. We show that Corollary 1 holds under the most unfavorable distribution and beta distributions.

4.2. Skewness versus semivariance

We let $q_{\mathrm{sv}}$ represent the candidate solution based on semivariance.

Corollary 2

In a special case with $m_1=0.5$ , $s=s_{\mathrm{k}}=0$ , and $c/p=0.5$ , the skew-based solution satisfies

$$q^{*}=\frac{(2m_2+\sqrt{4m_2-1})(\sqrt{4m_2-1}-4m_2+1)}{4m_2-1-(8m_2-3)\sqrt{4m_2-1}},$$

(21)

while the semivariance-based solution satisfies

$$q_{\mathrm{sv}}=0.5+0.25\sqrt{4m_2-1}>q^{*}\ge F^{-1}\left(\frac{p-c}{p}\right)=0.5$$

(22)

for any feasible distribution $F$ in set $\mathrm{\Omega }$ . Thus, $Z(q^{*}|F)$ dominates $Z(q_{\mathrm{sv}}|F)$ in this special case.

With the cost-to-price ratio being $0.5$ and this special case’s conditions on moments, the optimal inventory level is $q_F=0.5$ for any feasible distribution $F$ . With concavity, we know that if the inventory level is further away from $q_F=0.5$ , then the expected profit decreases. After verifying that $q^{*}$ is lower than $q_{\mathrm{sv}}$ , we conclude that the skew-based solution must outperform the semivariance-based solution in the special case. In Appendix C of the online E-Companion, we provide additional numerical examples to show that Corollary 2 can be generalized to the case with the cost-to-price ratio beyond $0.5$ .

5. Extensions

In this section, we first consider divergence-based solutions and then extend the analysis to a risk-averse newsvendor.

5.1. Divergence-based analysis

The divergence-based method has become increasingly popular (see Van Parys et al., 2021, for an updated literature review). This method typically starts with constructing an empirical distribution from the historical data. Let $\mathrm{\Theta }$ be the collection of realized demands according to the available history. We denote $P({\theta }_i)$ as the empirical probability mass for the realized ${\theta }_i$ .

5.1.1. KL divergence

We denote $Q({\theta }_i)$ as a feasible distribution satisfying the Kullback–Leibler (KL) divergence constraint

$$KL=\sum _{{\theta }_i\in \mathrm{\Theta }}P\left({\theta }_i\right)\ln \left(\frac{P\left({\theta }_i\right)}{Q\left({\theta }_i\right)}\right)\le \tau ,$$

(23)

where the parameter $\tau$ is chosen by the firm. We solve the robust optimization model

$$Z_{\tau }=\underset{q}{max}\underset{Q\left({\theta }_i\right)}{inf}\sum _{{\theta }_i\in \mathrm{\Theta }}Q\left({\theta }_i\right)Z({\theta }_i,q)$$

to determine the robust solution $q_{\tau }$ . The key difference between $Z$ and $Z_{\tau }$ is the ambiguity set. The former uses moments, but the latter uses the empirical distribution and KL divergence.

We observe that the extreme distribution associated with $q_{\tau }$ reassigns the probability mass by reducing the likelihood for higher realized demands but increasing that for lower ones. To be specific, nature optimizes the following Lagrangian to determine its best response $Q(\cdot )$ to any given $q$ :

$$\begin{array}{rl}\mathcal{L}=& -cq+p\sum _{{\theta }_i\in \mathrm{\Theta }}Q({\theta }_i)min({\theta }_i,q)\\ & +{\gamma }_1[\sum _{{\theta }_i\in \mathrm{\Theta }}P({\theta }_i)\ln \left(\frac{P\left({\theta }_i\right)}{Q\left({\theta }_i\right)}\right)-\tau ]\\ & +{\gamma }_2[\sum _{{\theta }_i\in \mathrm{\Theta }}Q({\theta }_i)-1].\end{array}$$

The FOC with respect to $Q({\theta }_i)$ satisfies

$$\begin{array}{rl}\frac{\mathrm{\partial }\mathcal{L}}{\mathrm{\partial }Q\left({\theta }_i\right)}& =pmin({\theta }_i,q)-{\gamma }_1\frac{P\left({\theta }_i\right)}{Q\left({\theta }_i\right)}+{\gamma }_2\\ & =0\text{for }\mathrm{\forall }{\theta }_i\in \mathrm{\Theta }.\end{array}$$

We recall that ${\gamma }_1\ge 0$ is the Lagrangian multiplier associated with the budget constraint (23), that ${\gamma }_2$ is the Lagrangian multiplier associated with the total probability constraint, and that ${\gamma }_2$ has an unrestricted sign. Knowing the empirical probability mass $P({\theta }_i)$ allows us to find the extreme distribution:

$$Q({\theta }_i)=\frac{{\gamma }_{1}P({\theta }_i)}{pmin({\theta }_i,q)+{\gamma }_2}.$$

(24)

Since the ex-post sales revenue $pmin({\theta }_i,q)$ is increasing in ${\theta }_i$ , equation (24) implies that $Q({\theta }_i)$ is decreasing in ${\theta }_i$ , ceteris paribus. In other words, nature must decrease (increase) the chance that higher (lower) values of realized demands can occur. This property of the extreme distribution explains why $q_{\tau }$ is more conservative than the moment-based solutions for larger $\tau$ because nature’s best response would reduce the mean (whereas in the moment-based model, nature keeps the mean unchanged).

Proposition 7

The optimal KL solution $q_{\tau }$ along with the two Lagrangian multipliers $({\gamma }_1,{\gamma }_2)$ can be characterized by the following three equations:

$$\begin{array}{rl}& \sum _{{\theta }_i\in \mathrm{\Theta }}\frac{{\gamma }_{1}P({\theta }_i)}{pmin({\theta }_i,q)+{\gamma }_2}=1\\ & \tau =\sum _{{\theta }_i\in \mathrm{\Theta }}P({\theta }_i)\ln \left(\frac{pmin({\theta }_i,q)+{\gamma }_2}{{\gamma }_1}\right)\end{array}$$

(25)

and

$$\frac{c}{p}=\sum _{{\theta }_i>q}Q({\theta }_i)=\sum _{{\theta }_i>q}\frac{{\gamma }_{1}P({\theta }_i)}{pmin({\theta }_i,q)+{\gamma }_2}.$$

(26)

Equations (25) and (26) are direct results of our joint optimization method, exemplifying its efficiency relative to the traditional two-stage approach.

5.1.2. Wasserstein distance

The Wasserstein distance is another popular divergence measure to construct the ambiguity set. Let $\sum P({\theta }_i)=\hat{P}({\theta }_i)$ be the cumulative distribution function of the empirical distribution $P(\cdot )$ . Similarly, let $\sum Q({\theta }_i)=\hat{Q}({\theta }_i)$ be the cumulative distribution function of the feasible distribution $Q({\theta }_i)$ . Let $N$ be the number of realized values in the empirical distribution. For convenience, we assume that ${\theta }_i<{\theta }_{i+1}$ , for $i=1,\dots ,N-1$ , meaning that the realized values ${\theta }_i$ are already ranked in increasing order. The Wasserstein distance equals

$$\begin{array}{rl}{W}_{\mathrm{d}}=& {\int }_{\mathrm{\Theta }}|\hat{P}(\theta )-\hat{Q}(\theta )|d\theta =({\theta }_2-{\theta }_1)|\hat{Q}({\theta }_1)-\hat{P}({\theta }_1)|\\ & +({\theta }_3-{\theta }_2)|\hat{Q}({\theta }_2)-\hat{P}({\theta }_2)|+\cdots \\ & +({\theta }_N-{\theta }_{N-1})|\hat{Q}({\theta }_{N-1})-\hat{P}({\theta }_{N-1})|.\end{array}$$

To determine the robust solution $q_{\mathrm{wd}}$ , we solve the robust optimization model

$$Z(q_{\mathrm{wd}})=\underset{q}{max}\underset{Q\left({\theta }_i\right)}{inf}\sum _{{\theta }_i\in \mathrm{\Theta }}Q\left({\theta }_i\right)Z({\theta }_i,q)$$

subject to the budget constraint $W_{\mathrm{d}}\le \tau$ and the total probability constraint $\sum _{{\theta }_i\in \mathrm{\Theta }}Q({\theta }_i)=1$ . We can find $q_{\mathrm{wd}}$ via the standard two-stage method, which involves a manageable computational effort, but finding our $q^{*}$ is almost effortless.

5.2. Conditional value at risk

The previous sections assume that the firm is risk-neutral. However, a firm may be risk-averse toward low profits. Therefore, in this subsection, we apply CVaR as the decision criterion to measure the firm’s level of risk-aversion. According to Rockafellar and Uryasev (2000), a risk-averse firm’s objective function under the CVaR criterion is

$$\begin{array}{rl}& CVa{R}_{\alpha }\left(Z(v,q,\theta )\right)\\ & =\underset{v,q}{max}\{v-\frac{1}{\alpha }E{(v-pmin(\theta ,q)+cq)}^{+}\},\end{array}$$

(27)

where $\alpha \in (0,1]$ is a parameter that measures risk preference and $v$ is the value at the risk level. If the demand distribution is known to be $F$ , then the $q$ that maximizes CVaR satisfies

$$q_{CR}=F^{-1}(\alpha \cdot \left(\frac{p-c}{p}\right)).$$

If $\alpha =1$ , then $q_{CR}$ is identical to the standard newsvendor solution. Thus, a smaller $\alpha$ leads to understocking.

Lemma 4

Given that the selling price $p$ ( $>c$ ) is fixed, if the firm is risk-averse and uses the CVaR criterion, then $v^{*}=(p-c)q.$

Hence, according to Lemma 4, we can rewrite equation (27) as

$$CVa{R}_{\alpha }\left(Z(q,\theta )\right)=\underset{q}{max}\{(p-c)q-\frac{1}{\alpha }{(pq-pmin(\theta ,q))}^{+}\}.$$

When the distribution $F$ is unknown, the firm’s objective function becomes

$$\begin{array}{c}Z(q,\theta )=\underset{q}{max}\underset{F}{inf}\{(p-c)q-\frac{1}{\alpha }{(pq-pmin(\theta ,q))}^{+}\}\\ \\ \text{s.t. }{\int }_0^{\mathrm{\infty }}{\theta }^{k}dF\left(\theta \right)=m_k,k=0,1,2,3.\end{array}$$

It is easy to verify that $Z(q,\theta )$ is continuous with respect to $\theta$ for any given $q$ and that $Z(q,\theta )$ is finite for any given $(q,\theta )$ . Thus, the dual model with known first two moments (mean and variance) is

$$\begin{array}{r}\begin{array}{l}Z=\underset{q,y_i}{max}\{y_0+y_1\mu +y_2({\mu }^2+{\sigma }^2)\}\\ \text{s.t. }{y}_0+y_1\theta +y_2{\theta }^2\le (p-c)q-\frac{1}{\alpha }{(pq-pmin(\theta ,q))}^{+},\mathrm{\forall }\theta .\end{array}\end{array}$$

(28)

Let $p^{\mathrm{\prime }}=(1/\alpha )p$ and $c^{\mathrm{\prime }}=c+((1-\alpha )/\alpha )p$ . The corresponding SIP model of equation (28) is

$$\begin{array}{l}Z=\underset{q,y_i}{max}\{y_0+y_1\mu +y_2({\mu }^2+{\sigma }^2)\}\\ \text{s.t. }\{\begin{array}{cc}{y}_0+y_1\theta +y_2{\theta }^2\le p^{\mathrm{\prime }}\theta -c^{\mathrm{\prime }}q,& \text{if}\theta \le q\\ y_0+y_1\theta +y_2{\theta }^2\le (p^{\mathrm{\prime }}-c^{\mathrm{\prime }})q,& \text{if}\theta >q,\end{array}\end{array}$$

which is similar to Scarf’s (1958) model (identical to Scarf’s (1958) model when $\alpha =1$ ). Thus, we have the following conclusions:

Corollary 3

When only the mean and standard deviation of the demand distribution are known, the firm’s most unfavorable distribution $F_{CVaR}$ satisfies

$$\begin{array}{l}Pr(\tilde{\theta }=m_1-\sqrt{\frac{p-\alpha (p-c)}{\alpha (p-c)}}\sqrt{m_2-m_1^2})=\frac{\alpha (p-c)}{p}\\ Pr(\tilde{\theta }=m_1+\sqrt{\frac{\alpha (p-c)}{p-\alpha (p-c)}}\sqrt{m_2-m_1^2})=1-\frac{\alpha (p-c)}{p}.\end{array}$$

The firm’s robust optimal order quantity is

$$q^{*}=m_1+\frac{1}{2}\sqrt{m_2-m_1^2}(\sqrt{\frac{\alpha (p-c)}{p-\alpha (p-c)}}-\sqrt{\frac{p-\alpha (p-c)}{\alpha (p-c)}}),$$

which is the mid-point of the two realized values. Additionally, the firm’s optimized utility equals $Z^{*}=(p-c)q^{*}$ .

Similarly, if the skewness is also known, then the SIP model is

$$\begin{array}{l}Z=\underset{v,q,y_i}{max}\{y_0+y_1\mu +y_2({\mu }^2+{\sigma }^2)+y_3{m}_3\}\\ \text{s.t.}\{\begin{array}{cc}{y}_0+y_1\theta +y_2{\theta }^2+y_3{\theta }^3\le p^{\mathrm{\prime }}\theta -c^{\mathrm{\prime }}q,& \theta \le q\\ y_0+y_1\theta +y_2{\theta }^2+y_3{\theta }^3\le (p^{\mathrm{\prime }}-c^{\mathrm{\prime }})q,& \theta >q\end{array}\end{array}$$

Applying the conclusions summarized in Propositions 1 to 6, we can obtain the firm’s most unfavorable distribution, the robust optimal order quantity, and the optimized utility, which are summarized in Appendix B of the E-Companion for reference.

6. Numerical study

In this section, we investigate the efficiency of the proposed robust strategy $q^{*}$ through numerical experiments on simulated data and actual industry data. Our focus is on whether skewness is a better measure than semivariance.

6.1. Simulated data

The characteristics of the true demand distribution affect the performance of any robust ordering policy. The beta distribution’s adaptability to a wide variety of scenarios and data types makes it highly useful in both simulations and practical applications. Readers can refer to the link https://en.wikipedia.org/wiki/Beta_distribution to see many different shapes of beta distributions (e.g., U-shaped, inverse U-shaped, positively or negatively skewed). Each beta distribution involves two parameters $(\alpha ,\beta )$ such that

$$\begin{array}{rl}{m}_1& =\frac{\alpha }{\alpha +\beta },m_2=\frac{\alpha (\alpha +1)}{(\alpha +\beta )(\alpha +\beta +1)},\text{and }\\ S_{\mathrm{k}}& =\frac{2(\beta -\alpha )\sqrt{\alpha +\beta +1}}{(\alpha +\beta +2)\sqrt{\alpha \beta }}.\end{array}$$

(29)

For each chosen beta distribution, we compute the semivariance by

$$\begin{array}{rl}SV& =\frac{\begin{array}{l}{\int }_{m_1}^1{(\theta -m_1)}^{2}f(\theta |\alpha ,\beta )d\theta -{\int }_0^{m_1}{(\theta -m_1)}^{2}f(\theta |\alpha ,\beta )d\theta \end{array}}{m_2-m_1^2},\end{array}$$

(30)

where $f(x|\alpha ,\beta )$ is the probability density function of the beta distribution. We find that $S_{\mathrm{k}}$ and $SV$ have the same sign for each pair $(\alpha ,\beta )$ . Additionally, for each given pair $(\alpha ,\beta )$ , we use the moment-related data in equations (29) and ( 30) to compute the information required to find $q^{*}$ and $q_{\mathrm{sv}}$ .

We define the distribution-dependent benchmark $q(\alpha ,\beta )$ by solving

$${\int }_0^{q(\alpha ,\beta )}f(x|\alpha ,\beta )dx=\frac{p-c}{p}$$

and then compute the benchmark profit $\overline{Z}(q(\alpha ,\beta ))$ as

$$\overline{Z}\left(q(\alpha ,\beta )\right)=-cq(\alpha ,\beta )+{\int }_0^{1}pmin(x,q(\alpha ,\beta ))f(x|\alpha ,\beta )dx.$$

To validate the impact of skewness on decision-making, we first select nine skewness values within the range of $[-2,+2]$ with a step length of $0.5$ . Next, we randomly draw 1,500 numbers between $0.1$ and $1$ as the $\alpha$ values. For each given skewness $S_{\mathrm{k}}$ , we use the $S_{\mathrm{k}}$ and $\alpha$ to calculate the corresponding $\beta$ values using

$$S_{\mathrm{k}}=\frac{2(\beta -\alpha )\sqrt{\alpha +\beta +1}}{(\alpha +\beta +2)\sqrt{\alpha \beta }}.$$

We limit $\alpha$ because if we expand the range of drawing $\alpha$ beyond $[0.1,1]$ , then the above equation may not yield any real solution for certain $S_{\mathrm{k}}$ . As we generate 1,500 $(\alpha ,\beta )$ pairs for each $S_{\mathrm{k}}$ value, we have a total of 13,500 combinations of $(\alpha ,\beta ,S_{\mathrm{k}})$ . We recall that the first four propositions show that the cost-to-price ratio plays a crucial role in determining the robust optimization solution. We select $p=5$ , $10$ , and $20$ , and set the production cost to $c=1$ so that the cost-to-price ratio takes the values of $\{0.2,0.1,0.05\}$ .

Table 2 devotes five rows to each of the three prices. The first two rows for each price list the percentage gap $\mathrm{\Delta }Z(q)=(1-Z(q)/\overline{Z}(q(\alpha ,\beta )))\times 100\mathrm{\%}$ for $q^{*}$ and for $q_{\mathrm{sv}}$ across all nine skewness values. A smaller gap indicates a better performance. To help us to determine if one policy dominates the other, the next three rows for each price count the number of $(\alpha ,\beta )$ pairs for which either candidate is better. Specifically, the rows denoted “>” list the number of instances for which $Z(q^{*}|\alpha ,\beta )>Z(q_{\mathrm{sv}}|\alpha ,\beta )$ , that is, the number (out of 1,500) for which $q^{*}$ strictly outperforms $q_{\mathrm{sv}}$ . The rows denoted “<” list the number of instances for which $Z(q^{*}|\alpha ,\beta )<Z(q_{\mathrm{sv}}|\alpha ,\beta )$ , that is, the number (out of 15,000) for which $q_{\mathrm{sv}}$ strictly outperforms $q^{*}$ .

Table 2.

Performance gaps of all candidate policies with simulated data.

$p$	Skewness	$-2$	$-1.5$	$-1$	$-0.5$	0	0.5	1	1.5	2	Avg.
4	$\mathrm{\Delta }Z(q^{*})$	0.00%	0.00%	0.06%	0.73%	0.93%	0.35%	2.37%	0.40%	0.70%	0.62%
	$\mathrm{\Delta }Z(q_{\mathrm{sv}})$	0.79%	1.30%	2.01%	2.48%	1.92%	0.51%	0.93%	8.03%	16.19%	3.80%
	>	1,500	1,500	1,500	1,284	1,267	570	0	1,500	1,500	1,180
	<	0	0	0	216	233	930	1,500	0	0	320
6	$\mathrm{\Delta }Z(q^{*})$	0.00%	0.00%	0.01%	0.16%	1.58%	0.54%	0.35%	1.99%	0.43%	0.56%
	$\mathrm{\Delta }Z(q_{\mathrm{sv}})$	0.32%	0.37%	0.65%	1.23%	1.81%	1.41%	0.38%	0.18%	1.73%	0.90%
	>	1,500	1,500	1,500	1,500	807	1,368	482	0	1,500	1,129
	<	0	0	0	0	693	132	1,018	1,500	0	371
8	$\mathrm{\Delta }Z(q^{*})$	0.00%	0.00%	0.00%	0.06%	0.72%	0.76%	0.34%	0.70%	2.37%	0.55%
	$\mathrm{\Delta }Z(q_{\mathrm{sv}})$	0.58%	0.59%	0.45%	0.56%	1.01%	1.31%	0.77%	0.13%	0.15%	0.62%
	>	1,500	1,500	1,474	1,498	1,014	1,217	1,362	192	0	1,084
	<	0	0	26	2	486	283	138	1,308	1,500	416

The bold values indicate the best results for each price.

We observe that $q^{*}$ outperforms $q_{\mathrm{sv}}$ , delivering a smaller average gap and a higher tally of “>” cases. Despite the better performance delivered by $q^{*}$ , we argue that such study (relying on underutilizing the available information) does not offer added value for practitioners because in practice the underlying demand distribution remains unknown. This motivates the study of the following subsection.

6.2. Industry data

We use data from the beverage company mentioned in section 1 to evaluate the performance of various robust production quantities. The company has 2 factories and 18 distribution centers, supplying 72 SKUs to external consumers. Since the numbers of available data points for each SKU are different, we choose $20$ SKUs that have similar numbers of observations without exhibiting a significant trend (so that the assumption of stationary demands holds). With production cost being scaled to $c=1$ , the prices vary so that the cost-to-price ratio $c/p$ takes three values from $\{0.2,0.1,0.05\}$ , consistent with subsection 6.1. For these discrete empirical distributions, we compute the semivariance as

$$SV=\frac{\begin{array}{l}\sum _{{\theta }_i>\overline{\theta }}{({\theta }_i-\overline{\theta })}^2-\sum _{{\theta }_i<\overline{\theta }}{({\theta }_i-\overline{\theta })}^2\end{array}}{\sum _{i=1}^n{({\theta }_i-\overline{\theta })}^2},$$

(31)

where $n$ is the number of days in the testing periods, ${\theta }_i$ is the real demand on day $i$ ( $i=1,2,\dots ,n$ ), $\overline{\theta }=(1/n)\sum _{i=1}^n{\theta }_i$ is the sample mean, and $(1/(n-1))\sum _{i=1}^n({\theta }_i-\overline{\theta }{)}^2$ is the sample variance. Although equation (31) is based on the observed sample data, other formulas for semivariance, including the unmarked equation on p. 3150 in Natarajan et al. (2018), are based on a known distribution $F$ .

Natarajan et al. (2018) used data from all the time periods to compute the moments and then implemented candidate policies in all periods. In contrast, we divide the data into training and testing periods. To be specific, we choose the first $T=30,60,90,120,150$ days as the training periods to estimate the moments that we use to construct various candidate policies. We then use the remaining days as the testing period to assess the performance of the proposed ordering policies.

This approach more closely matches a manager’s actual process. At whatever time the manager decides to estimate the demand’s distribution moments, he can use only the data known at that time, which he treats as training data. Certainly, the estimated moments based on these training data may differ from those of the testing data (i.e., the demand in the subsequent future months) in any practical implementation (especially when the underlying distribution is unknown). With these potentially inaccurate moments, our numerical study directly assesses the robustness (i.e., the reliable performance against ambiguous distributions) of various policies.

For any candidate inventory level $q$ , we compute the average daily profit as

$$Z\left(q\right)=\frac{1}{n}\sum _{i=1}^n[pmin({\theta }_i,q)-cq].$$

We test the performance of four policies: two moment-based solutions $q^{*}$ and $q_{\mathrm{sv}}$ , and two divergence-based solutions $q_{\tau }$ (KL divergence) and $q_{\mathrm{wd}}$ (Wasserstein distance). For $13$ different cost-to-price ratios and five different lengths of training period, we calculate the profit of $20$ selected SKUs under each policy during the test periods.

Next, we apply a retrospective optimization method based on full information to determine a constant inventory level that leads to the optimized retrospective profit $Z_{\mathrm{ro}}$ , which serves as the benchmark to evaluate the performance of each candidate policy. This retrospective method has been widely applied in the literature for a performance benchmark (e.g., Govindarajan et al., 2021; Natarajan et al., 2018). Specifically, we compute the optimized retrospective profit by solving

$$Z_{\mathrm{ro}}=\underset{q\ge 0}{max}\left\{\frac{1}{n}\sum _{i=1}^n[pmin({\theta }_i,q)-cq]\right\}.$$

(32)

Essentially, the retrospective method assumes that the demand history $\{{\theta }_i\}$ is known and then selects a constant inventory level to maximize the average daily profit. Hence, the retrospectively optimized profit $Z_{\mathrm{ro}}$ represents a loose upper bound. Nonetheless, we compute the percentage gap using

$$\mathrm{\Delta }Z(\cdot )=\frac{Z_{\mathrm{ro}}-Z(\cdot )}{Z_{\mathrm{ro}}}\times 100\mathrm{\%}.$$

(33)

Clearly, a smaller value of $\mathrm{\Delta }Z(\cdot )$ indicates a better performance of the candidate policy.

In addition to comparing $q^{*}$ and $q_{\mathrm{sv}}$ , we also apply the results of Section 5.1 with budgets of uncertainty $\tau \in \{0.01,0.05,0.1\}$ to construct divergence-based solutions $q_{\tau }$ and $q_{\mathrm{wd}}$ . However, the empirical distribution is also based on the training periods to maintain consistency with $q^{*}$ . To be specific, if the training periods include the first $30$ days, we use the demand data in the first $30$ days to estimate the moments and to establish the empirical distribution $P({\theta }_i)$ . Then, we construct both $q^{*}$ and $q_{\mathrm{sv}}$ using the closed-form results, but we implement a two-stage method to solve for $q_{\tau }$ and $q_{\mathrm{wd}}$ .

Table 3 reports the values of $\mathrm{\Delta }Z(\cdot )$ and uses bold numbers to indicate the smallest average performance gap among the four candidates for each combination of price and training period length. The fourth to eighth columns correspond to five different training period lengths, while the ninth column shows the average performance gap. The data in Table 3 suggests that the proposed $q^{*}$ policy is the best. Unsurprisingly, with longer training periods, the proposed $q^{*}$ policy performs better.

Table 3.

Performance gaps of various robust candidates.

			Training periods (%)
Price	Policy	Budget of uncertainty	30	60	90	120	150	Avg.
$p=4$	$\mathrm{\Delta }Z(q^{*})$	N/A	4.55	4.42	3.71	3.50	3.28	3.89
	$\mathrm{\Delta }Z(q_{\mathrm{sv}})$		6.78	7.57	6.15	6.05	8.68	7.05
	$\mathrm{\Delta }Z(q_{\tau })$	$\tau =0.01$	7.56	7.79	6.43	5.58	4.43	6.36
		$\tau =0.05$	14.30	6.05	9.62	6.60	6.72	8.66
		$\tau =0.10$	14.65	8.76	8.75	8.28	4.49	8.98
	$\mathrm{\Delta }Z(q_{\mathrm{wd}})$	$\tau =0.01$	5.54	5.07	3.89	3.20	2.68	4.08
		$\tau =0.05$	5.56	6.01	4.65	3.68	2.88	4.56
		$\tau =0.10$	5.81	5.80	5.35	4.67	4.61	5.25
$p=6$	$\mathrm{\Delta }Z(q^{*})$	N/A	3.77	3.09	3.01	2.24	1.71	2.76
	$\mathrm{\Delta }Z(q_{\mathrm{sv}})$		4.68	3.05	2.96	2.26	1.46	2.88
	$\mathrm{\Delta }Z(q_{\tau })$	$\tau =0.01$	9.70	8.02	6.46	4.61	4.04	6.57
		$\tau =0.05$	9.12	6.96	8.94	6.20	4.96	7.24
		$\tau =0.10$	12.78	7.35	9.62	7.32	4.01	8.21
	$\mathrm{\Delta }Z(q_{\mathrm{wd}})$	$\tau =0.01$	5.26	5.46	4.43	3.00	2.14	4.06
		$\tau =0.05$	5.73	5.42	4.48	3.17	2.43	4.25
		$\tau =0.10$	5.71	5.35	4.69	3.22	2.59	4.31
$p=8$	$\mathrm{\Delta }Z(q^{*})$	N/A	4.41	2.61	2.52	1.74	1.36	2.53
	$\mathrm{\Delta }Z(q_{\mathrm{sv}})$		5.15	2.86	2.91	2.15	1.34	2.88
	$\mathrm{\Delta }Z(q_{\tau })$	$\tau =0.01$	8.72	8.15	5.21	6.13	6.23	6.89
		$\tau =0.05$	9.00	7.62	6.56	6.98	6.66	7.36
		$\tau =0.10$	9.17	9.88	7.44	6.79	2.65	7.19
	$\mathrm{\Delta }Z(q_{\mathrm{wd}})$	$\tau =0.01$	4.94	4.12	4.14	3.11	2.24	3.71
		$\tau =0.05$	4.94	4.15	4.21	3.23	2.23	3.75
		$\tau =0.10$	5.02	4.12	4.35	3.33	2.23	3.81

To validate whether the out-performance is statistically significant, we also conducted a pairwise $t$ -test. Since all policies are implemented in the same set of SKU data, the pairwise $t$ -test is appropriate to detect any dominating performance. ⁴ We report the results of the pairwise $t$ -test in Table 4, in which we define the difference as the profit under our skewness policy minus that under the alternative policy. Thus, a positive $t$ -stat implies a positive difference (e.g., our skewness policy yields a higher profit).

Table 4.

Pairwise $t$ -test of industry data.

	Rival policy	$p$ Value	Paired difference $d$	$t$ -Stat	$df$	$SD$
$p=4$	$q_{\mathrm{sv}}$	$0.0005$	$68.8249$	$3.5946$	$99$	$191.4650$
	$KL$	$0.0007$	$347.8309$	$3.4875$	$99$	$997.3751$
	$WD$	$0.0062$	$97.2879$	$2.7967$	$99$	$347.8659$
$p=6$	$q_{\mathrm{sv}}$	$0.0003$	$156.5234$	$3.7824$	$99$	$413.8193$
	$KL$	$0.0000$	$702.0728$	$4.3920$	$99$	$\mathrm{1,598.5172}$
	$WD$	$0.0015$	$241.4715$	$3.2657$	$99$	$739.4184$
$p=8$	$q_{\mathrm{sv}}$	$0.0001$	$243.2397$	$4.0572$	$99$	$599.5215$
	$KL$	$0.0002$	$\mathrm{1,079.6224}$	$3.9335$	$99$	$\mathrm{2,744.7154}$
	$WD$	$0.0018$	$326.2077$	$3.2162$	$99$	$\mathrm{1,014.2513}$

Note. KL = Kullback–Leibler; WD = Wasserstein distance.

The analysis using industrial data is the closest to practice as we encounter several realistic challenges. First, the underlying distribution of the industrial data is unknown. We cannot preassign a distribution to evaluate the profit or from which to draw random numbers to form empirical estimates. Second, when learning from early periods, sampling errors exist, causing the estimated moments or the empirical distribution to differ from the true values. Third, the distribution drifts with unknown direction and magnitude. Despite these challenges, our skewness policy performs impressively.

6.3. Comparison

The crucial difference among various methods is the ambiguity set, which determines nature’s action space. First, moment-based ambiguity sets typically allow the realized demands to differ from those in the empirical distribution but assume the mean and variance to be constant. In the divergence-based ambiguity sets, the mean demand can either increase or decrease. Second, divergence-based methods use the empirical distribution as the starting point and then consider deviations in probability masses subject to budget constraints. Notably, nature varies only the probability masses without changing the support of the demand. For instance, if the empirical distribution does not include the point $\tilde{\theta }=7$ , then divergence-based ambiguity sets also exclude this point.

To illustrate the technical differences and consequences, we provide the following example based on a binomial distribution with parameters $(n,\pi )=(10,0.5)$ , so the input parameters include $m_1=5$ , $m_2=27.5$ , $m_3=162.5$ , and $S_{\mathrm{k}}=0$ . Let $p=2$ and $c=1$ , making the newsvendor ratio $(p-c)/p=0.5$ . Table 5 reports the results. Using Corollary 2 and scaling parameters, we confirm that the skewness policy yields $q^{*}=5.31$ and the semivariance policy yields $q_{\mathrm{sv}}=5.79$ . When using the two divergence-based methods, we set the budget of uncertainty to be $\tau =0.4$ and find that $q_{\tau }=q_{\mathrm{wd}}=4$ .

Table 5.

Performance of different policies without distribution drift (binomial demand).

Policy	$q_F$	$q^{*}$	$q_{\mathrm{sv}}$	$q_{\tau }$	$q_{\mathrm{wd}}$
Order quantity	$5$	$5.31$	$5.79$	$4$	$4$
Profit	$3.77$	$3.69$	$3.58$	$3.52$	$3.52$
$\frac{q-q_F}{q_F}\times 100\mathrm{\%}$	N/A	$6.2\mathrm{\%}$	$15.8\mathrm{\%}$	$-20\mathrm{\%}$	$-20\mathrm{\%}$
$\frac{Z(q_F|F)-Z(q|F)}{Z(q_F|F)}\times 100\mathrm{\%}$	N/A	$2.1\mathrm{\%}$	$5.0\mathrm{\%}$	$6.6\mathrm{\%}$	$6.6\mathrm{\%}$

As a benchmark, we consider no distribution drift (such that the empirical distribution perfectly predicts the underlying distribution $F$ ). Based on the given binomial distribution, the optimal solution is $q_F=5$ , yielding an optimized expected profit of $Z(q_F|F)=3.77$ . The performance gaps of the skewness policy, semivariance policy, and two divergence-based methods are respectively $2.1\mathrm{\%}$ , $5.0\mathrm{\%}$ , and $6.6\mathrm{\%}$ . Our skewness policy emerges as the best candidate. The ranking of the performance gaps is attributable to the gap in the order quantity. Relative to $q_F=5$ , the skewness policy overstocks by $6.2\mathrm{\%}$ , the semivariance rival overstocks by $15.8\mathrm{\%}$ , but the divergence-based methods understock by $20\mathrm{\%}$ .

To deepen our understanding of the causes of understocking, we investigate nature’s strategy in responding to $q_{\tau }$ and $q_{\mathrm{wd}}$ . We report the results in Table 6, in which the second column describes the empirical distribution, the third column characterizes nature’s best response to $q_{\tau }$ , and the fifth column identifies nature’s best response to $q_{\mathrm{wd}}$ .

Table 6.

Extreme distributions under $q_{\tau }$ and $q_{\mathrm{wd}}$ .

Demand	Empirical probability	Nature’s response to $q_{\tau }$	Divergence	Nature’s response to $q_{\mathrm{wd}}$	Wasserstein distance
$0$	$0.001$	$0.315$	$-0.006$	$0.201$	$0.200$
$1$	$0.010$	$0.024$	$-0.009$	$0.010$	$0.000$
$2$	$0.044$	$0.054$	$-0.009$	$0.044$	$0.000$
$3$	$0.117$	$0.096$	$0.023$	$0.117$	$0.000$
$4$	$0.205$	$0.126$	$0.099$	$0.015$	$0.190$
$5$	$0.246$	$0.152$	$0.119$	$0.246$	$0.000$
$6$	$0.205$	$0.126$	$0.099$	$0.205$	$0.000$
$7$	$0.117$	$0.072$	$0.057$	$0.117$	$0.000$
$8$	$0.044$	$0.027$	$0.021$	$0.044$	$0.000$
$9$	$0.010$	$0.006$	$0.006$	$0.000$	$0.010$
$10$	$0.001$	$0.001$	$0.000$	$0.000$	$0.000$
Total	$1.000$	$1.000$	$0.400$	$1.000$	$0.400$

Comparing the second and third columns of Table 6, we observe that nature’s behavior is consistent with equation (24). For example, nature drastically increases the original probability mass from $Pr(\tilde{\theta }=0)=0.001$ to $Pr(\tilde{\theta }=0)=0.315$ while decreasing the probability mass from $Pr(\tilde{\theta }=5)=0.246$ to $Pr(\tilde{\theta }=5)=0.152$ . As a result of diverting probability mass from high demands to low demands, nature’s strategy reduces the mean demand from $m_1=5$ to $m_1^{\mathrm{\prime }}=3.23$ . Similar observations can be made after comparing the second and fifth columns of Table 6. Under nature’s best response to $q_{\mathrm{wd}}$ , the mean demand decreases from $m_1=5$ to $m_1^{\mathrm{\prime }}=4.14$ .

We emphasize that the understocking problem in both $q_{\tau }$ and $q_{\mathrm{wd}}$ persists regardless of the empirical distribution or cost parameters. We explain this phenomenon from the zero-sum game perspective, in which nature wishes to minimize the firm’s profit by choosing a distribution from the ambiguity set. As we have mentioned, moment-based ambiguity sets preserve the value of the mean, but divergence-based ambiguity sets allow the mean to decrease. In Taylor’s expansion (C-1), we observe that the mean creates the first-order effect on the firm’s expected profit. Hence, nature (as the firm’s opponent) must try to reduce the mean demand, forcing the firm to lower the inventory level.

Next, we consider the impact of distribution drifts. First, suppose that the demand makes an upward drift to become a binomial distribution with $(n,\pi )=(10,0.53)$ (implying that the mean demand moderately increases by $6\mathrm{\%}$ from $5$ to $5.3$ ). We can verify that under the drifted distribution, the skewness policy $q^{*}=5.31$ yields an expected profit of $4.02$ ; the semivariance policy $q_{\mathrm{sv}}=5.79$ yields an expected profit of $3.97$ ; and both divergence-based policies $q_{\tau }=q_{\mathrm{wd}}=4$ yield an expected profit of $3.66$ . Thus, under an upward distribution drift, our skewness policy remains the best candidate. Finally, we suppose that the demand makes a downward drift. The drifted distribution is a binomial distribution with $(n,\pi )=(10,0.4)$ (implying that the mean demand drops significantly by $20\mathrm{\%}$ ). We can verify that under the drifted distribution, the skewness policy $q^{*}=5.31$ yields an expected profit of $2.32$ ; the semivariance policy $q_{\mathrm{sv}}=5.79$ yields an expected profit of $2.0$ ; and both divergence-based policies $q_{\tau }=q_{\mathrm{wd}}=4$ yield an expected profit of $2.8$ . In general, only when the underlying demand distribution makes a significant downward drift can divergence-based methods outperform our skewness policy.

Adjusting the budget of uncertainty $\tau$ is common among divergence-based methods since these methods implicitly allow the mean demand to decrease. A larger $\tau$ must expand nature’s action space by letting the mean further decrease. As a result, the firm must lower the inventory level in response to nature’s actions. For example, in a special case with $\tau =0$ , nature’s only option is the empirical distribution, making both $q_{\tau }$ and $q_{\mathrm{wd}}$ converge to the distribution-dependent optimal solution $q_F=5$ in the example shown in Table 6. When the budget increases from $0$ to $\tau =0.4$ , both $q_{\tau }$ and $q_{\mathrm{wd}}$ drop by $20\mathrm{\%}$ . Therefore, the adjustment of the budget $\tau$ causes a monotonic effect on the order quantities. Intuitively, if the distribution drift is believed to be weakly upward, then reduce $\tau$ ; if the distribution drift is believed to be downward, then increase $\tau$ . However, adjusting the budget level $\tau$ does not alter the fundamental trend that divergence-based methods understock.

6.4. Impact of inaccurate moments

To elaborate on the impact of inaccurate moments, we provide two numerical examples as follows. For expositional simplicity, we write $B(n,\pi )$ as the binomial distribution with probability mass $Pr(\tilde{\theta }=i)=(n!/(i!(n-i)!)){\pi }^i(1-\pi {)}^{n-i}$ . In the first example, we assume that the underlying distribution is $F=B(20,0.4)$ and the cost-to-price ratio is $c/p=0.6$ . The optimal order quantity is $q_F=7$ and the optimized profit is $Z_F=3.92$ . For each $n\in \{15,\dots ,25\}$ , we use the moments of $B(n,0.4)$ as the inputs to construct $q_n^{*}$ . Clearly, with $n<20$ ( $n>20$ ), the input data understates (overstates) all $m_i$ . After identifying $q_n^{*}$ , we evaluate the profit $Z(q_n^{*}|F)$ using the underlying distribution $F$ . In the second example, we assume that the underlying distribution is $F=B(20,0.6)$ and the cost-to-price ratio is $c/p=0.4$ . The optimal order quantity is $q_F=13$ and the optimized profit is $Z_F=15.88$ . We report the results in Figure 2(a) and (b), which display similar patterns.

Figure 2.

Impact of inaccurate moments. (a) Binomial distribution $(n,\pi )=(20,0.4)$ and $c/p=0.6$ and (b) binomial distribution $(n,\pi )=(20,0.6)$ and $c/p=0.4$ .

From Figure 2, we can conclude that when $n$ increases, the skewness order quantity $q_n^{*}$ increases as the input moments change from understatements to overstatements. However, the resultant profit $Z(q_n^{*}|F)$ is concave with respect to $n$ , with a peak that never exceeds $Z_F$ . The concave pattern suggests that when the accuracy of input data deteriorates, the performance of the skewness policy degrades.

7. Conclusion

The robust inventory management problem addresses stochastic inventory management for which the underlying demand distribution is unknown. Motivated by practical needs and theoretical curiosity to go beyond the first two moments, we develop a robust inventory model based on the first three moments (i.e., mean, variance, and skewness). We demonstrate that information about skewness should be ignored by the firm only when the demand is unbounded and the cost-to-price ratio is high. Therefore, we recommend that firms acquire information about skewness if either (i) the maximum demand can be reliably estimated or (ii) the cost-to-price ratio is low. If neither of these two conditions holds, then firms can use the mean and variance (without skewness) in their planning process. For the unbounded case, we specify the threshold of the cost-to-price ratio that determines whether using skewness leads to more profitable decisions. We also find that if using skewness is beneficial under the most unfavorable distribution, then the firm’s optimized utility increases with respect to variance but decreases with respect to skewness. This result contradicts the economics and finance literature. Numerical studies (based on real data with unknown distributions from a large beverage company) also show that our moment-based policy, on average, delivers a more consistent performance than one based on the Kullback–Leibler divergence, one based on Wasserstein distance, and one based on semivariance. Given the convenience of estimating skewness (i.e., using Excel), we recommend that firms include skewness in their robust inventory management (Table 7).

Table 7.

When does skewness matter?

Ratio $(c/p)$ demand	Bounded	Unbounded
Low	Skewness matters	Skewness matters
High	Skewness matters	Variance matters

Methodologically, we contribute by developing a joint optimization procedure to derive closed-form optimal solutions. This method is more efficient than traditional methods, which require solving the second-stage problem but often encounter intractable FOCs. Our analytical results find the firm’s most unfavorable distribution—which is essential for robust optimization—in four jointly exhaustive cases for bounded demands. For unbounded demand, these results require two cases that depend on whether using skewness provides greater profit.

Our joint optimization method paves the way for several other applications in which practitioners recognize the importance of including higher moments. The remarkable advantage of the joint optimization method is that we bypass intermediate roadblocks and directly find the final answer. Thus, our joint optimization approach is a parsimonious process for solving these complex and important problems.

Supplemental Material

sj-pdf-1-pao-10.1177_10591478261446418 - Supplemental material for When does skewness matter in robust inventory management?

Supplemental material, sj-pdf-1-pao-10.1177_10591478261446418 for When does skewness matter in robust inventory management? by Feng Tao, Yao-Yu Wang, Zhaolin Li and H Neil Geismar in Production and Operations Management