The News Vendor Model

All code in one file

Introduction

Here are some python implementations of the newsvendor model as discussed in [hopp08:_factor_physic]. We illustrate the code with some examples. We use the following notation:

\[\begin{split}\begin{split} X &= \text{The stochastic demand}, \\ G &= \text{The distribution of } X,\\ Q &= \text{The ordering quantity}, \\ c_s &= \text{shortage cost}, \\ c_o &= \text{overage cost}. \end{split}\end{split}\]

Computation of Expected Profit

Here we show how to compute the expected profit for various examples of the newsvendor model. To see that the expected profit must be defined as

\[Z(Q) = (c_s + c_o) E\min\{Q,X\} - c_o Q,\]

where \(E(Y)\) represents the expected value of the random variable \(Y\), we reason as follows. Suppose the newsvendor does not sell anything, which occurs when \(X=0\). This must, clearly, cost \(c_o Q\), since his entire ordered quantity \(Q\) will yield a return value of \(c_o\) per unit returned. Supppose now that the newsvendor sells some amount \(X\leq Q\). Each unit sold will yield a revenue equal to the shortage cost \(c_s\). Moreover, each unit sold will not be returned, hence result in an extra yield \(c_o\). Finally, since the newsvendor can never sell more than \(Q\) items, the amount sold should be capped with taking the minimum of \(Q\) and \(X\).

import numpy as np
import scipy.stats as stats
from scipy.stats import poisson, norm

c_s = 10.
c_o = 4.

Suppose demand is Poisson distributed with mean \(\mu\). The following generates the probability mass function.

mu = 10
G = poisson(mu)
p = G.pmf(np.arange(100))

The expected profit is now a straightforward implementation of the formula for \(Z(Q)\). Use that the expected number of sold items can be written as

\[ES = E \min\{Q, X\} = \sum_{i=0}^\infty \min\{Q, i\} p_i,\]

where \(p_i = P(X=i)\), i.e., the probability mass function at \(i\).

def Z(Q):
    ES = sum((min(Q,i)*p[i] for i in range(len(p))))
    return (c_s+c_o)*ES - c_o*Q

Thus, we implement the above definition of \(Z(Q)\).

It can be shown that optimal order quantity is such that \(G(Q^*) = c_s/(c_s+c_o)\). Note that we have to round this result to an integer.

Qstar = int(np.ceil(G.ppf(c_s/(c_s+c_o))))

print Qstar
12

Lets see whether this is correct

print "Q, profit"
for Q in range(Qstar-3, Qstar+5):
    print Q, Z(Q)
Q, profit
9 78.8956110016
10 82.484594999
11 84.3220384963
12 84.5671724481
13 83.4853817786
14 81.3828798619
15 78.551298483
16 75.2336641292

Indeed, at \(Q = Q^*\) the maximum profit is attained.

Value of Perfect Information

Suppose it would be possible to predict precisely the demand. Then the profit of the newsvendor would equal \(c_s E X\). Thus, the, so-called, value of perfect information is

\[c_s E X - Z(Q),\]

where we assume that \(Q\) is taken to be the optimal order quantity.

def P():
    expectedDemand = sum((i*p[i] for i in range(len(p))))
    return c_s*expectedDemand

print P() -  Z(Q), P()/Z(Q)
24.7663358708 1.32919220614

This is the expected profit that can be made by taking the effort to obtain perfect information about demand. Interestingly, this is not that much, provided the newsvendor uses the optimal order quantity.

A Christmas Lights Example

This is another example of the newsvendor problem. Here we implement the newsvendor algorithms in a class.

def expect(X, f, lb, ub):
    if hasattr(X, 'dist'):
        return X.dist.expect(f, lb = lb, ub = ub, *X.args, **X.kwds)
    else:
        return X.expect(f, lb = lb, ub = ub)

For some technical reasons we need this function. For the techically inclined, currently rv_frozen objects in scipy.stats do not support the expect method. This function provides a work around for this deficiency.

from scipy.integrate import quad
import pylab as pl

class newsVendor(object):
    def __init__(self, co, cs, X):
        self.co = co
        self.cs = cs
        self.X = X # demand distribution during lead time
        self.G = self.X.cdf

    def unitsOver(self, Q):
        return expect(self.X, lambda x: Q-x, 0, Q)

    def unitsShort(self, Q):
        return expect(self.X, lambda x: x-Q, Q, np.inf)

    def Y(self,Q):
        return self.co*self.unitsOver(Q) + self.cs*self.unitsShort(Q)

    def Qstar(self):
        return self.X.ppf(self.cs/(self.cs+self.co))

    def profit(self, Q):
        investedValue = self.co * Q
        expectedSales = expect(self.X, lambda x: np.minimum(x,Q), 0, np.inf)
        return (self.co + self.cs)*expectedSales - investedValue

Christmas lights

X = norm(loc = 10000., scale = 1000.)
news = newsVendor(co = 2.5, cs = 5., X= X)
Q =  news.Qstar()

print Q
print news.Y(Q)
print news.Y(10440) # rounding inaccuracy in FP
10430.7272993
2726.99831006
2727.11539102

Selling cakes

Yet another example. Suppose you would try to sell cakes at the entrance of some university building. The basic problem is to estimate how much profit you expect to make. In other words, is it, business-wise, an interesting case to try to sell cakes like this?

To obtain some quantitative understanding of this problem we make some simplyfing assumptions. First, cakes cannot be carried over to the next day, hence the scrap value is zero. The buying price from the local backery \(p_b = 25\) cents. Suppose you would sell the cakes for \(p_s = 1\) Euro. Under this price assumption we make the following guess about the distribution of the sales \(X\):

\[\begin{split}\begin{align*} P(X \leq 50) &= 2/10 \\ P(50 < X \leq 100) &= 3/10 \\ P(100 < X \leq 150) &= 3/10 \\ P(150 < X \leq 200) &= 3/20 \\ P(200 < X \leq 300) &= 1/20 \\ \end{align*}\end{split}\]

Since we don’t have any further information, we just assume that within each interval the demand is uniformly distributed. So now the question is: what is the optimal buying quantity \(Q^\star\) from the backery?

First we need the overage and shortage cost per cake.

buy = 0.25
sell = 1.
scrap = 0.
co = buy - scrap
cs = sell - buy

Now we build the demand distribution.

thresholds = [50, 100, 150, 200, 300]
p=[2,3,3,1.5,0.5]

probs = []
grid = []
t0 = 0
for t,q in zip(thresholds, p):
     probs += [q]*(t-t0)
     grid += range(t0, t)
     t0 = t

tot = sum(probs)
probs = [p/tot for p in probs]

X =  stats.rv_discrete(values=(grid,probs), name='demand')

The optimal order quantity:

news = newsVendor(co = co, cs = cs, X= X)
Q =  news.Qstar()
print Q
print news.profit(Q)
147.0
59.7566666667

So we expect to earn roughly 60 euro a day with this business.

Let’s plot the expected profit as a function of \(Q\).

Z = []
Qvalues = range(300)
for Q in Qvalues:
    Z.append(news.profit(Q))
pl.plot(Qvalues,Z)
../../_images/newsvendor_figure13.png

How sensitive is the expected profit on the distribution? To get some insight into this, let’s replace the above demand by normally distributed demand with the same mean and standard deviation.

mean, std = X.mean(), X.std()
print mean, std
X = stats.norm(loc = mean, scale = std)
news = newsVendor(co = co, cs = cs, X= X)
Q =  news.Qstar()
print Q
print news.profit(Q)
110.214285714 67.2902510856
155.60087036
62.706926496

Finally, suppose it would be possible to obtain complete information about the demand. In such cases we would know precisely what we would sell the next day. Hence, this would be the amount of cakes to order the day before. In such a case our earnings would be \(c_s E (X)\), i.e., the expected demand times the shortage cost.

The value of complete information is then the difference between what we can maximally make, that is, \(c_o E (X)\), and the earnings we make when we would use the optimal newsvendor policy, that is, the policy with the optimal order quantity \(Q^\star\).

print cs*X.mean()
print news.profit(Q)
82.6607142857
62.706926496

We can give this difference an interesting interpration. It is clear that it costs money/effort to obtain better predictions of the demand. When this cost exceeds 20 Euro on a daily basis, it simply does not make sense to invest in obtaining it.