Problem 3.5.17

Walther has a roadside produce stand where he sells oats, peas, beans, and barley. He buys these products at wholesale prices of, respectively, $1.05, $3.17, $1.99, and $0.95; he sells them at per-pound retail prices of, respectively, $1.29, $3.76, $2.23, and $1.65. Each day the amount demanded (in pounds) could be as little as zero for each product, and as much as 10, 8, 14, and 11 for oats, peas, beans, and barley; respectively; he sells only whole-pound amounts, no partial pounds. Assume a discrete uniform distribution for daily demand for each product over its range; assume as well that Walther always has enough inventory to satisfy all demand. The summer selling season is 90 days. Create a spreadsheet simulation that will, for each day as well as for the whole season, simulate Walther’s total cost, total revenue, and total profit.

n <- 90
day <- matrix(NA, n, 7)

for (i in 1:n) {
  day[i, 1] <- floor((10 - 0) * runif(1)) # oats
  day[i, 2] <- floor((8 - 0) * runif(1)) # peas
  day[i, 3] <- floor((14 - 0) * runif(1)) # beans
  day[i, 4] <- floor((11 - 0) * runif(1)) # barley
  day[i, 5] <- sum(day[i, 1:4] * c(1.05, 3.17, 1.99, 0.95))
  day[i, 6] <- sum(day[i, 1:4] * c(1.29, 3.76, 2.23, 1.65))
  day[i, 7] <- day[i, 6] - day[i, 5]
  }

t(data.frame(mean_cost = mean(day[, 5]), 
             mean_revenue = mean(day[, 6]), 
             mean_profit = mean(day[, 7]),
             total_cost = sum(day[, 5]), 
             total_revenue = sum(day[, 6]), 
             total_profit = sum(day[, 7])))

##                     [,1]
## mean_cost       34.37422
## mean_revenue    42.63622
## mean_profit      8.26200
## total_cost    3093.68000
## total_revenue 3837.26000
## total_profit   743.58000

import numpy as np
import pandas as pd

n = 90
day = np.empty(shape=[n, 7])

for i in range(0, n):
  day[i, 0] = int(np.random.uniform(0, 10, 1)) # oats
  day[i, 1] = int(np.random.uniform(0, 8, 1)) # peas
  day[i, 2] = int(np.random.uniform(0, 14, 1)) # beans
  day[i, 3] = int(np.random.uniform(0, 11, 1)) # barley
  day[i, 4] = sum(day[i, 0:4] * [1.05, 3.17, 1.99, 0.95])
  day[i, 5] = sum(day[i, 0:4] * [1.29, 3.76, 2.23, 1.65])
  day[i, 6] = day[i, 5] - day[i, 4]

df = pd.DataFrame(np.array([np.mean(day[:, 4]), np.mean(day[:, 5]), np.mean(day[:, 6]), 
  np.sum(day[:, 4]), np.sum(day[:, 5]), np.sum(day[:, 6])]), [['mean_cost', 
  'mean_revenue', 'mean_profit', 'total_cost', 'total_revenue', 'total_profit']])

print df

##                          0
## mean_cost        34.032556
## mean_revenue     42.400444
## mean_profit       8.367889
## total_cost     3062.930000
## total_revenue  3816.040000
## total_profit    753.110000

Simio Simulation

This example is a static simulation. Simio is an “special-purpose dynamic simulation software.” (p. 39). Therefore, although it can probably be hacked, a Simio simulations are not designed for this type of scenario.

Problem 3.5.1

Extend the simulation of throwing two dice in Section 3.2.1 in each of the following ways (one at a time, not cumulatively):

1. Instead of 50 throws, extend the spreadsheet to make 500 throws, and compare your results. Tap the F9 key to get a “fresh” set of random numbers and thus a fresh set of results.

2. Load the dice by changing the probabilities of the faces to be something other than uniform at \(1/6\) for each face. Be care to ensure that the probabilities sum to 1.

3. Use @RISK (see Section 3.2.3), or another Excel add-on for static Monte Carlo spreadsheet simulation, to make 10,000 throws of the pair of fair dice, and compare your results to the true probabilities of getting the sum equal to each of 2, 3, …, 12 as well as to the true expected value of 7.

sum_table <- matrix(1:6, nrow=6, ncol=6)
sum_totals <- sum_table + t(sum_table)
round(ftable(c(sum_totals)) / length(sum_totals), 3)

##      2     3     4     5     6     7     8     9    10    11    12
##                                                                   
##  0.028 0.056 0.083 0.111 0.139 0.167 0.139 0.111 0.083 0.056 0.028

After the Excel add-in @RISK was installed, the formula =RiskUniform(1,6) was entered in cells A1 and B1, and the formula =RiskOutput("@RISK Sum") + SUM(A1:B1) was enetered in cell C1. Then, under the “@RISK” tab, 10,000 simulations were specified and “Start Simulation” was clicked to generate the results below. The mean is exactly 7 and the probability of each sum value in the histogram appears to be in the range of the exact probabilities.

Problem 3.5.5

In the Monte Carlo integration of Section 3.2.2, add to the spreadsheet calculation of the standard deviations of the 50 individual values, and use that, together with the mean already in cell H4, to compute the 95% confidence interval on the exact integral in cell I4; does your confidence interval contain, or “cover” the exact integral? Repeat all this but with a 90% confidence interval and then with a 99% confidence interval.

\[\bar { x } \pm { z }_{ { \alpha }/{ 2 } }\frac { s }{ \sqrt { n } }\]

The confidence intervals were calculated in Excel using AVERAGE(\({ \left\{ (b-a)\phi _{ \mu ,\sigma }(x_{ 1 }) \right\} }_{ i=1 }^{ n }\)) \(\pm\) NORM.INV(\(\alpha/2\), \(\mu\), \(\sigma\)) * STDEV.S(\({ \left\{ (b-a)\phi _{ \mu ,\sigma }(x_{ 1 }) \right\} }_{ i=1 }^{ n }\)) / SQRT(\(n\)). The 95% and 90% confidence intervals cover the exact integral. The 99% confidence interval does not.

Problem 3.5.17

See above.