# The News Vendor Model¶

## 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:

## 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

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

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

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\):

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)
```

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.