Monte Carlo Simulation | Blockchain Consultants

Part 6: Comparing Alternative Scenarios

In the first and second articles of this series, we defined a Monte Carlo Simulation (MCS) as “a sampling experiment whose aim is to estimate the distribution of a quantity of interest that depends on one or more stochastic input variables”. It is a quantitative technique whose modern development was based on the work of scientists’ John von Neumann, Stanislav Ulam, and Nicholas Metropolis in the Manhattan Project for the development of the first atomic bomb during World War II. Conceptually, it is a method based on statistics for determining the likelihood of multiple outcomes based on repeated random sampling.

Simulation techniques in general, and MCS in particular, excel in comparing alternative scenarios. The ability to model a lot of assumptions about the system under study makes them the most powerful and flexible operation research tools.

One of the most important purposes of any simulation study is to compare system alternatives using quantitative statistical measures which describe the performance of the system. The final goals are to determine which scenario performs the best according to the decision problem and to gain insight into the system’s behavior.

So, an appropriate simulation methodology must include a thorough selection of inputs, objective and measurable outcome metrics, and a final analysis and ranking of alternative scenarios.

We will illustrate these ideas and concepts for a Monte Carlo Simulation through a classical example in operation management models: the news vendor inventory problem.

The news vendor inventory problem appears frequently in business and operation research decision problems. It is modeled by the so-called Newsvendor Model. It tries to answer the following question: how much to order for a single product that is perishable and characterized by a random demand and a fixed selling price. The vendor wants to maximize his/her profit but, as the product must be available before the selling period, he/she must decide the size of the order before observing the real product demand.

As detailed in the first article of this series, the general setting of the problem is as follows:

• There is a vendor, a supplier, and customers.
• The demand is uncertain for an upcoming period of time (a single-period model).
• There is a decision variable (the inventory level Q) that must be computed to determine the one that maximizes the profit (the optimal decision).
• The uncertain customer demand is represented by a random variable d with its corresponding probability density function f(d) and cumulative distribution function F(d).
• The seller perfectly knows the fixed initial cost per unit (c), the selling price per unit (s), and the salvage value (u). The salvage value is the estimated value of a product at the end of its useful life. The habitual business situation indicates that the selling price is higher than the purchase cost, and this in turn is higher than the salvage value.
• Each unit of demand above the inventory level is lost as a sale; each unit of demand below the inventory level is sold at the salvage value. If Q > d, (Q — d) units are left over and must be salvaged by the vendor; otherwise, (d — Q) units represent lost sales.

Several extensions to the classical newsvendor problem have been developed. One of them involves advertising-sensitive demand: how customer demand is influenced by advertising and marketing efforts.

This extension models the customer demand as a function of the marketing efforts. A certain number of potential customers receive the advertisement, but the number who will make a purchase decision will depend on a factor named advertising effectiveness. On the other hand, there is the so-called advertising intensity related to the number of potential customers to whom the marketing effort will be directed. This parameter will determine the advertising and marketing costs according to the formula below:

where α is a constant.

So, our newsvendor model with advertising-sensitive demand includes the following equations:

1) Cost_per_Order = c * Q

2) Revenue_per_Sales = s * minimum(d, Q)

3) Revenue_per_Salvage = u * maximum(0, Q — d)

4) Advertising Cost = C(AI)

and the profit equation is:

Profit = Revenue_per_Sales + Revenue_per_Salvage — Cost_per_Order — Advertising Cost

So, the profit is a function of the order quantity (inventory level) and the customer demand (advertising intensity). We’ll use a Monte Carlo Simulation to estimate the order quantity that maximizes the expected value of the profit for different values of the advertising intensity. This calculation will be performed by sampling over given data and scenarios. Every combination of an order quantity and an advertising intensity is a scenario.

We are going to simulate a hypothetical newsvendor problem using data adopted from Negahban [1].

Our company would like to know how many products to order from the supplier for selling them in the upcoming season.

The supplier offers us four alternatives:

· Alt1: 500 units at a cost of $200/unit

· Alt2: 1000 units at a cost of $190/unit

· Alt3: 1500 units at a cost of $175/unit

· Alt4: 2000 units at a cost of $165/unit

We defined, for each alternative, the corresponding selling price and assumed, initially, a salvage value equivalent to 50% of the forecasted selling value.

Finally, we chose six different advertising intensities (5%, 10%, 15%, 20%, 25%, and 30%). The following table shows, for each one, the corresponding values of the loc and scale parameters (mean and standard deviation) for the normal distribution that represents the random demand for the different intensities.

We have four supplier options and six advertising intensities: this translates into 24 scenarios. We are going to replicate 3000 runs for each scenario, so the half-width of the confidence interval for our principal point estimator (profit) should be reasonable for our decision problem.

The python code for our hypothetical news vendor problem is as follows:

@author: darwt# Import Modules
import pandas as pd
import numpy  as npfrom scipy import stats
from scipy.stats import semimport matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
from prettytable import PrettyTable
your_path="your_path"#..................................................................
# initialization module
list_Price_Per_Unit   = [330, 300, 270, 230]
list_Cost_Per_Unit    = [200, 190, 175, 165]multiplier = 0.50
list_Salvage_Per_Unit = list(np.array(list_Price_Per_Unit) * multiplier)list_advert_intens = [5, 10, 15, 20, 25, 30]
list_Mean_Normal = [288, 938, 1306, 1691, 1940, 2243]
list_Std_Normal  = [211, 262, 281,  308,  294,  289]
length1 = len(list_advert_intens)list_of_orders = [500,1000,1500,2000]
length2 = len(list_of_orders)alfa = 0.5 
confidence = 0.95                      ## selected by the analyst 
# .......................................................
column_labels = ["Advert.Intens.", "Order Quantity",
"Mean","Std.Dev","Var.","Std. Error", 
"Median", "Skewness", "Kurtosis",                    
"CI Half Width", "CI LL", 'CI UL']df = pd.DataFrame(columns=column_labels)

The logic behind the MCS for this particular newsvendor model is described in the following lines of code:

listM1, listM2, listM3, listM4, listM5, listM6 =[],[],[], [], [], [] 
list_of_Means = [listM1, listM2, listM3, listM4, listM5, listM6]Number_of_Replications = 3000for i in range(length1):
for j in range(length2):
list_of_profits = []
for run in range(Number_of_Replications):
Qty_Demanded = np.random.normal(                                          loc=list_Mean_Normal[i],scale=list_Std_Normal[i],size=1) 
if Qty_Demanded >= 0:
Qty_Demanded = int(Qty_Demanded)
else:
Qty_Demanded = 0
Qty_Ordered  = list_of_orders[j]
Qty_Sold = np.minimum(Qty_Demanded, Qty_Ordered)
Qty_Left = np.maximum(0, Qty_Ordered - Qty_Demanded)
Revenue_per_Sales = Qty_Sold * list_Price_Per_Unit[j]            
Revenue_per_Salvage= Qty_Left * list_Salvage_Per_Unit[j]
Cost_per_Order =  Qty_Ordered * list_Cost_Per_Unit[j]
Cost_per_Advert = alfa * np.log(-1/((list_advert_intens[i]/100)-1))            
Cost_per_Advert = Cost_per_Advert * 100000
Profit = Revenue_per_Sales + Revenue_per_Salvage - Cost_per_Order - Cost_per_Advert
list_of_profits.append(Profit)

We used Numpy and Scipy to obtain a set of key descriptive statistical measures, particularly the expected value of the profits:

        media = np.mean(list_of_profits)
stand = np.std(list_of_profits)
var   = np.var(list_of_profits) 
std_error = sem(list_of_profits)
median = np.median(list_of_profits)
skew   = stats.skew(list_of_profits)
kurt   = stats.kurtosis(list_of_profits)
dof  = Number_of_Replications - 1    
t_crit = np.abs(stats.t.ppf((1-confidence)/2,dof))
half_width =  round(stand *t_crit/np.sqrt(Number_of_Replications),2)  
inf = media - half_width
sup = media + half_width  
inf = round(float(inf),2)
sup = round(float(sup),2)
list_of_statistics = []
list_of_statistics.append(list_advert_intens[i])
list_of_statistics.append(Qty_Ordered)
list_of_statistics.append(round(media,2))
list_of_statistics.append(round(stand,2))
list_of_statistics.append(round(var,2))
list_of_statistics.append(round(float(std_error),2))
list_of_statistics.append(round(median,2))
list_of_statistics.append(round(float(skew),2))
list_of_statistics.append(round(float(kurt),2))
list_of_statistics.append(round(half_width,2))
list_of_statistics.append(round(inf,2))
list_of_statistics.append(round(sup,2))
df.loc[len(df)] = list_of_statistics
list_of_Means[i].append(round(media,2))

We used Matplotlib to show a table with six rows and four columns for the 24-scenarios profits.

row_list = [str(x) + ' %'     for x in list_advert_intens]
col_list = [str(x) + ' Units' for x in list_of_orders]
subtitle_text="Salvage Value equal to %s" %multiplier +
' of the price value'fig, ax = plt.subplots(1,1)
ax.axis('tight')
ax.axis('off')
profits_table = ax.table(cellText  =  list_of_Means, 
rowLabels = row_list,
colLabels = col_list,                      
rowColours =["skyblue"]*(length1),  
colColours =["cyan"]*length2, 
cellLoc="center", loc="center",
bbox = [0.1, 0, 1.9, 1.0])
ax.set_title(" Profits according to Order Quantity & Advertising Intensity",fontsize=18, y= 1.2 , pad = 4)plt.figtext(0.5, 0.95,
subtitle_text,
horizontalalignment="center",
size= 12, style="italic",
color="black"
)profits_table.auto_set_font_size(False)
profits_table.set_fontsize(12)
for (row, col), cell in profits_table.get_celld().items():
if (row == 0) or (col == -1):
cell.set_text_props(
fontproperties=FontProperties(weight="bold"))plt.savefig(your_path +'Scenarios.png',bbox_inches="tight", dpi=150)
plt.show()

Then, we used Matplotlib for plotting the same data for better visualization:

fig, ax = plt.subplots()
for i in range(len(list_of_Means[0])):
plt.plot(list_advert_intens,[pt[i] for pt in list_of_Means],
label = list_of_orders[i],
linestyle="--", marker="*")
plt.legend()
plt.axhline(y = 0.0, color="k", linestyle="-")
fig.suptitle('Monte Carlo Simulation',  fontsize=20)
plt.xlabel('Advertising Intensity (%)', fontsize=16)
plt.ylabel('Profit (U$S)', fontsize=16)
plt.xticks()plt.savefig(your_path +'ScenarioChart.png',
bbox_inches="tight", dpi=150)
plt.show()

Finally, we used PrettyTable for the statistic report corresponding to the scenario with the maximum profit.

max_profit = df.loc[df["Mean"].idxmax()]
t = PrettyTable(['Statistic', 'Value'])
t.add_row(['Runs', Number_of_Replications])
t.add_row(['Advert.Intens.',max_profit['Advert.Intens.']])
t.add_row(['Order Quantity',max_profit['Order Quantity']])
t.add_row(['Mean', max_profit['Mean']])
t.add_row(['Median', max_profit['Median']])
t.add_row(['Variance', max_profit['Var.']])
t.add_row(['Stand. Error', max_profit['Std. Error']])
t.add_row(['Skewness', max_profit['Skewness']])
t.add_row(['Kurtosis', max_profit['Kurtosis']])
t.add_row(['Half Width',max_profit['CI Half Width']])
t.add_row(['CI inf', max_profit['CI LL']])
t.add_row(['CI sup', max_profit['CI UL']])
print(t)

Analysis

We made 3000 replications for each of the 24 scenarios, appending the statistical measures to a dataframe. Table 1 shows the mean profit for each scenario. The best numeric decision involves purchasing 1500 units and performing a marketing effort with a 25% of advertising intensity.

Figure 1 shows the same data in linear plots. The chart shows that we can reduce the advertising intensity to 20% without seriously harming our profits.

Table 2 shows the half-width of the confidence interval for the maximum profit. The half-width of the confidence interval habitually indicates the precision of the estimate of the measure of performance in a simulation study. The half-width shown in the table meets our precision requirements.

Finally, we repeated the simulation with two different percentages for the salvage value: 25% and 75 % of the initially forecasted selling value. Table 3 shows the mean profit for the new scenarios (75%). It can be seen that there are no changes in our decision problem: we must purchase 1500 units and make marketing efforts with a 20% advertising intensity.

This is another practical implementation of a Monte Carlo Simulation to solve a complex decision problem in a very simple manner. Undoubtedly, it is a very powerful technique to model, analyze, and compute performance metrics for many real-world systems.

Don’t forget to give a TIP, particularly if you add the article to a list.

[1] Negahban, A. “A hybrid simulation framework for the newsvendor problem with advertising and viral marketing”. Proceedings of the 2013 Winter Simulation Conference. R. Pasupathy, S.-H. Kim, A. Tolk, R. Hill, and M. E. Kuhl, eds.

Monte Carlo Simulation Republished from Source https://towardsdatascience.com/monte-carlo-simulation-bf31bb78d39c?source=rss—-7f60cf5620c9—4 via https://towardsdatascience.com/feed

<!–

–>