C. Simulation, part 2

Mon, Oct 24, 2022 Demand estimation is a technique to model how consumer behaviour changes due to changes in the price of the product, consumer income, or any other variable that impacts demand.1 The previous post demonstrated simple simulation through iterating through a normal distribution.

Demand estimation is a technique to model how consumer behaviour changes due to changes in the price of the product, consumer income, or any other variable that impacts demand.¹

The previous post demonstrated simple simulation through iterating through a normal distribution. It also went over what a confidence interval is for, and how it can be reduced by increasing sample size.

C.2 Demand Estimation

To use demand estimation for optimization (remember optimization?), the steps are:

Fit a distribution (normal, Weibull, Logistic, etc.) and its parameters such as mean and standard deviation to actual historical data. If there is no historical data, use a normal distribution and take an educated guess as to the mean and standard deviation.
Predict the optimal future outcome on actual input parameters. For example, the optimal number of widgets to sell based on the current variable cost of production and the current selling price.
Simulate over several iterations, changing the input parameters. Record the optimal solution. For example, change the variable costs and current selling prices for each iteration. Make note of the optimal input parameters.

C.2.1 Example with the Newsvendor Problem

Let’s put this simple 3-step process to work.

In the case of the 140-year-old Newsvendor Problem², let us assume over the last 30 days the daily number of newspapers sold at a busy location were as such: $$ \scriptsize{\color{#717171}{ num\_sold = [92, 96, 75, 90, 67, 88, 64, 78, 93, 78, 42, 115, 94, 81,\\ 114, 110, 92, 83, 107, 63, 68, 114, 93, 129, 97, 77, 61, 42, 81, 114] }} $$

Fit the distribution. In Excel, use the @Risk Fit function. In Python, StackOverflow user Saulo G.P. Castro demonstrated the equivalent³, modified by Timothy Davenport, and simplified below:

#!/usr/bin/env python

import warnings

import numpy as np
import scipy.stats as st
from scipy.stats._continuous_distns import _distn_names

num_sold=[
    92, 96, 75, 90, 67, 88, 64, 78, 93, 78, 42, 115, 94,
    81, 114, 110, 92, 83, 107, 63, 68, 114, 93, 129, 97,
    77, 61, 42, 81, 114
]
BINS = 200  # for the histogram used in the fit
SKIP_LIST = ['studentized_range']

best_distributions = []

# Get histogram of original data
y, x = np.histogram(num_sold, bins=BINS, density=True)
x = (x + np.roll(x, -1))[:-1] / 2.0

for dist_name in _distn_names:
    if dist_name in SKIP_LIST:
        continue
    try:
        with warnings.catch_warnings():
            warnings.filterwarnings('ignore')
            distribution = getattr(st, dist_name)
            params = distribution.fit(num_sold)

            # separate parts of parameters
            arg = params[:-2] # estimates for shape parameters
            loc = params[-2]  #  mean
            scale = params[-1]  # standard deviation

            # probability distribution function
            pdf = distribution.pdf(x, *arg, loc=loc, scale=scale)

            # SSE is the sum of squares for error. Really
            # important machine learning. The lowest SSE is
            # the best result.
            sse = np.sum(np.power(y - pdf, 2.0))

    except ValueError:
        pass
    best_distributions.append(
        (distribution.name, loc, scale, sse)
    )

best_distributions = sorted(
    best_distributions, key=lambda x:x[3]  # The lowest SSE wins.
)
best_distributions = [ round(i, 4) for i in best_distributions]

print(best_distributions[0])

And we get:

./best_dist.py
('dgamma', 9.3808, 85.7010, 0.23712)

The “dgamma” is the most reasonable distribution, mean 86, and standard deviation 9.38. If one explores the other results’ SSE’s, however, one will see that the competing distributions are all very close to each other. Using “norm” is not only viable, it is simpler. I will continue to use “dgamma” just for my own curiosity. Since it is close to a normal distribution, I will set its alpha (not required for norm), to 1.

Optimize with actual input variables. In this case, we want to optimize the vendor profits. We want specifically to know how many newspapers to buy from the publisher (stock) based on current publisher price and current retail price. Let us say the vendor buys at $1.00/pc, sells at $1.50/pc, and returns unsold papers at $0.40/pc. Since the mean of sales is 86 papers, we’ll generate stock values from 50 to 130. We will use the distribution to generate the uncertainty of actual purchases.

Note on code: @Risk requires the number of simulations to correspond to the range of decision variables. It’s just the way spreadsheets work. Here we can be more sophisticated since Python is general programming, but we won’t. Note also the SEED. This is used on pseudorandom numbers so random results are actually reproducible by third parties.

#!/usr/bin/env python

import numpy as np
import scipy.stats as st

# The decision variable is the stock quantities
START = 50
STOP = 130
stock_quantities = range(START, STOP+1, 5)  # Every 5.
buy_price = 1.00
sell_price = 1.50
sellback_price = 0.40

# The input variable is random, generated from the previous fit
sigma = 9.3808
mu = 86

# The output random variable is profit. We'll record all here.
profits = []

# Seed
SEED = 42  #  Life, the universe, and everything

# Simulate!
NUM_SIMULATIONS = int((STOP-START)/5)  # 16

# Generate a series of purchases based on distribution
random_purchases = st.dgamma.rvs(
    a=1, loc=mu, scale=sigma, random_state=SEED,
    size=int(NUM_SIMULATIONS)
)  # alpha=1 means no right skew, i.e. centered

i = 0
while i < NUM_SIMULATIONS:
    num_bought = stock_quantities[i]

    # You cannot sell more than you bought
    num_sold = min(num_bought, int(random_purchases[i]))
    num_unsold = num_sold - num_bought

    cost = buy_price*num_bought
    revenue = sell_price*num_sold + sellback_price*num_unsold
    profits.append((num_bought, round(revenue - cost, 2)))
    i += 1

print(profits)
sorted_profits = sorted(profits, key=lambda x:x[1])
print("Best result:", sorted_profits[-1])

The results:

[(50, 25.0), (55, 27.5), (60, 30.0), (65, 32.5), (70, 35.0),
(75, 37.5), (80, 40.0), (85, 42.5), (90, 45.0), (95, 47.5),
(100, 17.7), (105, 27.8), (110, 24.6), (115, 0.5),
(120, -21.7), (125, -15.4)]
Best result: (95, 47.5)

Our vendor should buy 95 newspapers. Notice, however, that that’s at the edge of a cliff. Profits go down rapidly. One might mess around with the seed. I tried seed values of 36 and 236. Ordering 90 and 95 newspapers seems to be optimal in each case, but ordering 90 is safer, as sometimes 95 is a steep drop in profits. Going over 105 risks negative profit but sometimes it is very profitable indeed – perhaps this is too volatile for our newsvendor!

Change the input parameters. What happens if the buy price goes to $1.75 and the sell price gets jacked to $2.25? Assume losing 20% of the readership, and the buyback only goes to $0.50 because the publisher is squeezing hard.

One will need to find the mean and standard deviation again since the number of sales will have changed.

...
buy_price = 1.75
sell_price = 2.25
sellback_price = 0.50

# The input variable is random, generated from the previous fit

# Original values
num_sold=[
    92, 96, 75, 90, 67, 88, 64, 78, 93, 78, 42, 115, 94,
    81, 114, 110, 92, 83, 107, 63, 68, 114, 93, 129, 97,
    77, 61, 42, 81, 114
]
num_sold_reduced = [int(0.8 * i) for i in num_sold]
mu = int(np.mean(num_sold_reduced))
sigma = np.std(num_sold_reduced)

...

And we get:

[(50, 25.0), (55, 27.5), (60, 30.0), (65, 32.5), (70, -14.5),
(75, 10.0), (80, -9.5), (85, 15.0), (90, 12.0), (95, 42.0),
(100, -49.0), (105, -16.25), (110, -19.25), (115, -74.5),
(120, -127.0), (125, -105.25)]
Best result: (95, 42.0)

Once again, ninety-five newspapers, but that’s a weird bump among all the losses. Changing the seed seemed to give the same answer: a low number of stock makes a consistent but low profit, and only a very small sweet spot gives a better profit, surrounded by loss.

One had better hope the publisher does not treat the newsvendors like this.

C.2.2 Experiment and conclusion

Mess around with the data. Although the example used Python’s double gamma distribution, the numbers given were very close to a normal distribution as well (which is why I set alpha to 1). The code samples are good for exploration and learning the concepts, but they are certainly no match for the ease of @Risk on Excel. A nice programming project would be to encapsulate the code in classes and methods, and add matplotlib calls to make pretty graphs.