Kelton et al. (2014). Simio and Simulation: Modeling, Analysis, Applications (3rd Ed.)

3.5.17) Walther has a roadside produce stand where he sells oats, peas, beans, and barley. He buys these products at per-pound 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.33, 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, and demand each day is independent of demand on other 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.

library(ggplot2)
library(dplyr)
library(tidyr)

produce <- matrix(NA, nrow = 90, ncol = 4)
colnames(produce) <- c("oats", "peas", "beans", "barley")

crp <- as.data.frame(matrix(NA, nrow = 90, ncol = 3))
colnames(crp) <- c("cost", "revenue", "profit")

cost <- c(1.05, 3.17, 1.99, 0.95)
revenue <- c(1.29, 3.76, 2.33, 1.65)
n_max <- c(10, 8, 14, 11)

sim <- as.data.frame(matrix(NA, nrow = 1000, ncol = 7))
colnames(sim) <- c(colnames(produce), colnames(crp))

for (n in 1:nrow(sim)) {
  for (i in 1:4) {
    produce[, i] <- sample(0:n_max[i], 90, replace = TRUE, 
                         prob = rep((1 / length(0:n_max[i])), length(0:n_max[i])))
    sim[n, i] <- sum(produce[, i])
  }
  crp$cost <- as.numeric(produce %*% cost)
  crp$revenue <- as.numeric(produce %*% revenue)
  crp$profit <- crp$revenue - crp$cost
  
  sim[n, ]$cost <- sum(crp$cost)   
  sim[n, ]$revenue <- sum(crp$revenue)
  sim[n, ]$profit <- sum(crp$profit)
}

knitr::kable(produce[1:10, ], caption = "Produce demand for first 10 days of 1000th simulated summer")
Produce demand for first 10 days of 1000th simulated summer
oats peas beans barley
7 8 0 5
8 7 9 10
0 8 9 2
6 2 0 10
0 5 0 6
7 7 4 5
1 2 6 10
8 3 6 5
10 0 3 2
4 8 11 7 
knitr::kable(sim %>% summarize(Avg_Cost = mean(cost),
                               Avg_Revenue = mean(revenue),
                               Avg_Profit = mean(profit)),
             caption = "Average values for 1000 simulated summer seasons")
Average values for 1000 simulated summer seasons
Avg_Cost Avg_Revenue Avg_Profit
3339.106 4220.069 880.9633 
res <- sim %>% 
  mutate(summer = row_number()) %>% 
  select(summer, cost:profit) %>%
  gather(name, value, cost:profit)

res$name <- factor(res$name, levels = c("cost", "revenue", "profit"))

ggplot(res) +
  geom_histogram(aes(value, fill = name), alpha = 0.6) +
  scale_fill_manual(name = "",
    values = c("cost" = "red", "revenue" = "blue", "profit" = "green"),
    labels = c("Cost", "Profit", "Revenue")) +
  labs(x = "Dollars", 
       y = "Count", 
       title = paste("Distribution of Total Cost, Revenue & Profit across",
                     "1000 Simulated 90-day Summer Seasons",
                     sep = "\n")
  )

ggplot(res) +
  geom_line(aes(summer, value, color = name, group = name), alpha = 0.6) +
  scale_color_manual(name = "",
    values = c("cost" = "red", "revenue" = "blue", "profit" = "green"),
    labels = c("Cost", "Profit", "Revenue")) +
  labs(x = "Summer", 
       y = "Dollars", 
       title = "Total Cost, Revenue & Profit for each Simulated Summer Season")