2  The newsvendor problem

Author

Affiliation

Aditya Mahajan

McGill University

Keywords

stochastic optimization, newsvendor problem

Image credit: https://americangallery.wordpress.com/category/cafferty-james-h/

Summary

The newsvendor problem is a simple model of stochastic optimization problem where a decision has to be made when there is uncertainty about the outcome. It also shows that for some stochastic optimization problems it is possible to obtain the qualitative properties of the nature of optimal solution.

Each morning, a newsvendor has to decide how many newspapers to buy before knowing the demand during the day. The newsvendor purchases a newspaper at a cost of $\mathrm{\$}p$ per newspaper and sells them at a cost of $\mathrm{\$}q$ per newspaper, where $q>p$. Any unsold newspapers at the end of the day have no salvage value.

Let $a$ denote the number of newspapers bought and $W$ denotes the demand. If $W<a$, then the newsvendor will sell $W$ newspapers and receive a total earnings of $qW-pa$. If $W\ge a$, then the newsvendor will sell $a$ newspapers and receive a total earning of $qa-pa$. Thus, the reward is $r(a,W)$, where

$$r(a,w)=\{\begin{array}{ll}qw-pa,& \text{if }w<a,\\ qa-pa,& \text{if }w\ge a.\end{array}$$

2.1 Interlude with continuous version

The problem above has discrete action and discrete demand. To build intuition, we first consider the case where both the actions and demand are continuous. Let $f(w)$ denote the probability density of the demand and $F(w)$ denote the cumulative probability density. Then, the expected reward is $$\begin{array}{}\text{(1)}& J(a)={\int }_0^a[qw-pa]f(w)dw+{\int }_a^{\mathrm{\infty }}[qa-pa]f(w)dw.\end{array}$$

To fix ideas, we consider an example where $p=0.5$, $q=1$, and the demand is a Kumaraswamy distribution with parameters $(2,5)$ and support $[0,100]$. The performance of a function of action is shown below.

p = 0.5
q = 1
r = function(w,a){ if(w<=a) { return q*w - p*a } else { return q*a - p*a } }

a_opt = inverseCDF( (q-p)/q )

config = ({
  // Kumaraswamy Distribution: https://en.wikipedia.org/wiki/Kumaraswamy_distribution
  a: 2,
  b: 5,
  max: 100
})

pdf = {
  const a = config.a
  const b = config.b

  return function(x) {
    var normalized = x/config.max
    return a*b*normalized**(a-1)*(1 - normalized**a)**(b-1)
  }
}

inverseCDF= {
  const a = config.a
  const b = config.b
  return function(y) {
     // Closed form expression for inverse CDF of Kumaraswamy distribution 
     return config.max * (1 - (1-y)**(1/b))**(1/a)
  }
}

points = { 
  const n = 1000
  var points = new Array(n)
  var cdf = 0;
  for (var i = 0 ; i < n; i++) {
    var x = config.max * i/n
    cdf = cdf + pdf(x) * 1/n
    points[i] = {demand: x, pdf: pdf(x), CDF: cdf, reward: r(x,action) }
  }
  return points
}

cost_values = {
  const n = 1000
  var points = new Array(n)
  for (var i = 0 ; i < n; i++) {
    var x = config.max * i/n
    points[i] = {action: x, performance: J(x)}
  }
  return points
}

J = {
  const a = config.a
  const b = config.b

  return function(action) {
    const n = 1000
    var cost = 0
    var w = 0
    for (var i = 0; i < n; i++) {
      w = config.max*i/n
      if (w <= action) {
        cost += (q*w - p*action)*pdf(w)/n
      } else {
        cost += (q*action - p*action)*pdf(w)/n
      }
    }
    return cost
  }
}

p = 0.5

q = 1

r = ƒ(w, a)

a_opt = 35.979082354039534

config = Object {a: 2, b: 5, max: 100}

pdf = ƒ(x)

inverseCDF = ƒ(y)

points = Array(1000) [Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, …]

cost_values = Array(1000) [Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, …]

J = ƒ(action)

viewof action = Inputs.range([0, 100], {value: 45, step: 0.01, label: "a"})

cost = Math.round(J(action)*100)/100

a

action = 45

cost = 10.51

plotJ = Plot.plot({
    grid: true,
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.dot([ [action, J(action)] ], {fill: "blue", r:4}),
    Plot.dot([ [a_opt, J(a_opt)] ], {fill: "red", r:4}),
    Plot.line(cost_values, {x:"action", y:"performance"})
  ]
})

plotPDF = Plot.plot({
  grid: true,
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.line([ [action,0], [action, pdf(action)] ], {stroke: "blue"}),
    Plot.line(points,{x:"demand", y:"pdf"}),
    Plot.areaY(points.filter(pt => pt.demand <= action),{x:"demand", y:"pdf", fill: "lightblue"}),
    Plot.areaY(points.filter(pt => pt.demand > action),{x:"demand", y:"pdf", fill: "pink"})
  ]
})

plotReward = Plot.plot({
  grid: true,
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.line(points, {x:"demand", y:"reward"})
  ]
})

(a) Performance: 10.51

(b) PDF of demand

(c) reward (as a function of demand)

Figure 2.1: An example to illustrate the results. Plot (a) shows the performance as a function of action; the blue dot shows the value of chosen action and the red dot shows the value of optimal action. Plot (b) shows the PDF of the demand where the blue shaded region shows the probability of getting a demand less than ordered goods and the red shaded region shows the probability of getting a demand greater than ordered goods. Plot (c) shows the reward function $r(a,\cdot )$, which depends on the values of $p$ and $q$.

In Figure 2.1(a), the plot of $J(a)$ is concave. We can verify that this is true in general.

Verify that $J(a)$ is concave

To verify that the function $J(a)$ is concave, we compute the second derivative: $$\frac{d^{2}J(a)}{d{a}^2}=-pf(a)-(q-p)f(a)=-qf(a)\le 0.$$

This suggests that we can use calculus to find the optimal value. In particular, to find the optimal action, we need to compute the $a$ such that $dJ(a)/da=0$.

Proposition 2.1 For the newsvendor problem with continuous demand, the optimal action is $$a=F^{-1}(1-\frac{p}{q}).$$ In the literature, the quantity $1-(p/q)$ is called the critical fractile.

Proof

Leibniz integral rule

$${\displaystyle \frac{d}{dx}}({\int }_{p(x)}^{q(x)}f(x,t)dt)=f(x,q(x))\cdot {\displaystyle \frac{d}{dx}}q(x)-f(x,p(x))\cdot {\displaystyle \frac{d}{dx}}p(x)+{\int }_{p(x)}^{q(x)}{\displaystyle \frac{\partial }{\partial x}}f(x,t)dt.$$

Using the Leibniz integral rule, the derivative of the first term of $\text{(1)}$ is $$[qa-pa]f(a)+{\int }_0^a[-p]f(w)dw=[qa-pa]f(a)-pF(a).$$

Similarly, the derivative of the second term of $\text{(1)}$ is $$-[qa-pa]f(a)+{\int }_a^{\mathrm{\infty }}(q-p)f(w)dw=-[qa-pa]f(a)+(q-p)[1-F(a)].$$

Combining the two, we get that $${\displaystyle \frac{dJ(a)}{da}}=-pF(a)+(q-p)[1-F(a)].$$

Equating this to $0$, we get $$F(a)={\displaystyle \frac{q-p}{q}}\text{or}a=F^{-1}(1-{\displaystyle \frac{p}{q}}).$$

Graphical interpretation of the result

The result of Proposition 2.1 has a nice graphical interpretation. Draw the CDF of the demand. The optimal action is the point where the CDF intersects the horizontal line $1-p/q$.

viewof pcts = Inputs.range([0.01,1], {value: 0.5, step: 0.01, label: "p"})
viewof qcts = Inputs.range([0.01,1], {value: 1, step: 0.01, label: "q"})

p

pcts = 0.5

q

qcts = 1

F_opt = 1 - pcts/qcts
plotOptCts = Plot.plot({
  grid: true,
  y: {domain: [0,1]},
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.line(points,{x:"demand", y:"CDF"}),
    Plot.line([ [0,F_opt], [points[points.length-1].demand, F_opt] ], {stroke: "red"})
  ]
})

F_opt = 0.5

2.2 Back to discrete version

Now, we come back to the problem with discrete actions and discrete demand. Suppose $W$ takes the values $\mathcal{W}=\{w_1,w_2,\dots ,w_k\}$ (where $w_1<w_2<\cdots <w_k$) with probabilities $\{{\mu }_1,{\mu }_2,\dots ,{\mu }_k\}$. It is easy to see that in this case the action $a$ should be in the set $\{w_1,w_2,\dots ,w_k\}$.

To fix ideas, we repeat the above numerical example when $\mathcal{W}=\{0,1,\dots ,100\}$.

pointsD = { 
  const n = 100
  var points = new Array(n)
  var cdf = 0
  for (var i = 0 ; i < n; i++) {
    var x = config.max * i/n
    cdf = cdf + pdf(x) * 1/n
    points[i] = {demand: x, pdf: pdf(x) * 1/n, CDF: cdf, reward: r(x,actionD) }
  }
  return points
}

cost_valuesD = {
  const n = 100
  var points = new Array(n)
  for (var i = 0 ; i < n; i++) {
    var x = config.max * i/n
    points[i] = {action: x, performance: J(x)}
  }
  return points
}


a_optD = Math.round(a_opt*100)/100

pointsD = Array(100) [Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, …]

cost_valuesD = Array(100) [Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, Object, …]

a_optD = 35.98

viewof actionD = Inputs.range([0, 100], {value: 45, step: 0.01, label: "a"})

costD = Math.round(J(actionD)*100)/100

a

actionD = 45

costD = 10.51

plotJD = Plot.plot({
    grid: true,
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.dot([ [actionD, J(actionD)] ], {fill: "blue", r:4}),
    Plot.dot([ [a_optD, J(a_optD)] ], {fill: "red", r:4}),
    Plot.line(cost_valuesD, {x:"action", y:"performance", curve: "step-after"})
  ]
})

plotPDFD = Plot.plot({
  grid: true,
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    // Plot.line([ [actionD,0], [actionD, pdf(actionD)] ], {stroke: "blue"}),
    Plot.line(pointsD,{x:"demand", y:"pdf", curve:"step-after"}),
    Plot.areaY(pointsD.filter(pt => pt.demand <= actionD),{x:"demand", y:"pdf", curve: "step-after", fill: "lightblue"}),
    Plot.areaY(pointsD.filter(pt => pt.demand >= actionD),{x:"demand", y:"pdf", curve: "step-after", fill: "pink"})
  ]
})

plotRewardD = Plot.plot({
  grid: true,
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.line(pointsD, {x:"demand", y:"reward", curve: "step-after"})
  ]
})

(a) Performance: 10.51

(b) PDF of demand

(c) reward (as a function of demand)

Figure 2.2: An example to illustrate the results. Plot (a) shows the performance as a function of action; the blue dot shows the value of chosen action and the red dot shows the value of optimal action. Plot (b) shows the PDF of the demand where the blue shaded region shows the probability of getting a demand less than ordered goods and the red shaded region shows the probability of getting a demand greater than ordered goods. Plot (c) shows the reward function $r(a,\cdot )$, which depends on the values of $p$ and $q$.

In the discrete case, the brute force search is easier (because there are a finite rather than continuous number of values). We cannot directly use the ideas from calculus because functions over discrete domain are not differentiable. But we can use a very similar idea. Instead of checking if $dJ(a)/da=0$, we check the sign of $J(w_{i+1})-J(w_i)$.

Proposition 2.2 Let $\{M_i{\}}_{i\ge 1}$ denote the cumulative mass function of the demand. Then, the optimal action is the largest value of $w_i$ such that $$M_i\le 1-\frac{p}{q}.$$

Proof

The expected reward for choice $w_i$ is $$\begin{array}{rl}J(w_i)& =\sum _{j<i}{\mu }_j[q{w}_j-p{w}_i]+\sum _{j\ge i}{\mu }_j[q{w}_i-p{w}_i]\\ & =-p{w}_i+q[\sum _{j<i}{\mu }_j{w}_j+\sum _{j\ge i}{\mu }_j{w}_i].\end{array}$$

Thus, $$\begin{array}{rl}J(w_{i+1})-J(w_i)& =-p{w}_{i+1}+q[\sum _{j<i+1}{\mu }_j{w}_j+\sum _{j\ge i+1}{\mu }_j{w}_{i+1}]\\ & +p{w}_i-q[\sum _{j<i}{\mu }_j{w}_j+\sum _{j\ge i}{\mu }_j{w}_i]\\ & =-p(w_{i+1}-w_i)+q[\sum _{j\ge i+1}{\mu }_j(w_{i+1}-w_i)]\\ & =(-p+q[1-M_i])(w_{i+1}-w_i).\end{array}$$ Note that $$M_i\le {\displaystyle \frac{q-p}{q}}⟺-p+q[1-M_i]\ge 0.$$ Thus, for all $i$ such that $M_i\le (q-p)/q$, we have $J(w_{i+1})\ge J(w_i)$. On the other hand, for all $i$ such that $M_i>(q-p)/q)$, we have $J(w_{i+1})<J(w_i)$. Thus, the optimal amount to order is the largest $w_i$ such that $M_i\le (q-p)/q$.

Graphical interpretation of the result

The structure of the optimal solution is the same for continuous and discrete demand distributions. In particular, the result of Proposition 2.2 has the same graphical interpretation as that of Proposition 2.1:

Draw the CDF of the demand. The optimal action is the point where the CDF intersects the horizontal line $1-p/q$.

viewof pdis = Inputs.range([0.01,1], {value: 0.5, step: 0.01, label: "p"})
viewof qdis = Inputs.range([0.01,1], {value: 1, step: 0.01, label: "q"})

p

pdis = 0.5

q

qdis = 1

F_opt_dis = 1 - pdis/qdis
plotOptDts = Plot.plot({
  grid: true,
  y: {domain: [0,1]},
  marks: [
    // Axes
    Plot.ruleX([0]),
    Plot.ruleY([0]),
    // Data
    Plot.line(pointsD,{x:"demand", y:"CDF", curve:"step-after"}),
    Plot.line([ [0,F_opt_dis], [pointsD[pointsD.length-1].demand, F_opt_dis] ], {stroke: "red"})
  ]
})

F_opt_dis = 0.5

Exercises

Exercise 2.1 (Qualitative properties of optimal solution) Intuitively, we expect that if the purchase price of the newspaper increases but the selling price remains the same, then the newsvendor should buy less newspapers. Formally prove this statement.

Hint: The CDF of a distribution is a weakly increasing function.

Exercise 2.2 (Monotonicity of optimal action) Consider two scenarios for the case with continuous demand and actions. In scenario 1, the demand is distributed according to PDF $f_1$. In scenario 2, it is distributed according to PDF $f_2$. Suppose $F_1(w)\le F_2(w)$ for all $w$. Show that the optimal action $a_1$ for scenario 1 is greater than the optimal action $a_2$ for scenario 2.

Hint: Plot the two CDFs and try to interpret the optimal decision rule graphically.

Exercise 2.3 (Selling random wind) The amount $W$ of power generated by the wind turbine is a positive real-valued random variable with probability density function $f$. The operator of the wind turbine has to commit to provide a certain amount of power in the day-ahead market. The price of power is $\mathrm{\$}p$ per MW.

If the operator commits to provide $a$ MWs of power and the wind generation $W$ is less than $a$, then he has to buy the balance $a-W$ from a reserves market at the cost of $\mathrm{\$}q$ per unit, where $q>p$. Thus, the reward of the operator is $r(a,W)$ where $$r(a,w)=\{\begin{array}{ll}pa,& \text{if }w>a\\ pa-q(a-w),& \text{if }w<a.\end{array}$$

Find the value of commitment $a$ that maximizes the expected reward.

Notes

Perhaps the earliest model of the newsvendor problem appeared in Edgeworth (1888) in the context of a bank setting the level of cash reserves to cover demands from its customers. The solution to the basic model presented above and some of its variants was provided in Morse and Kimball (1951); Arrow et al. (1952); Whitin (1953). See Porteus (2008) for an accessible introduction.

The property $F_1(w)\le F_2(w)$ used in Exercise 2.2 is called stochastic dominance. Later in the course, we will study how stochastic dominance is useful to establish monotonicity properties of general MDPs.

The example of selling random wind in Exercise 2.3 is taken from Bitar et al. (2012).

Arrow, K.J., Harris, T., and Marschak, J. 1952. Optimal inventory policy. Econometrica 20, 1, 250–272. DOI: 10.2307/1907830.

Bitar, E., Poolla, K., Khargonekar, P., Rajagopal, R., Varaiya, P., and Wu, F. 2012. Selling random wind. Hawaii international conference on system sciences, IEEE, 1931–1937.

Edgeworth, F.Y. 1888. The mathematical theory of banking. Journal of the Royal Statistical Society 51, 1, 113–127. Available at: https://www.jstor.org/stable/2979084.

Morse, P. and Kimball, G. 1951. Methods of operations research. Technology Press of MIT.

Porteus, E.L. 2008. Building intuition: Insights from basic operations management models and principles. In: D. Chhajed and T.J. Lowe, eds., Springer, 115–134. DOI: 10.1007/978-0-387-73699-0.

Whitin, S. 1953. The theory of inventory management. Princeton University Press.