Lab 08
CME Group Foundation Business Analytics Lab
Mona Tariq
4/3/2019
About
In this lab we will focus on sensitivity analysis and Monte Carlo simulations.
Sensitivity analysis is the study of how the uncertainty in the output of a mathematical model or system (numerical or otherwise) can be apportioned to different sources of uncertainty in its inputs. We will use the lpSolveAPI R-package as we did in the previous lab.
Monte Carlo Simulations utilize repeated random sampling from a given universe or population to derive certain results. This type of simulation is known as a probabilistic simulation, as opposed to a deterministic simulation.
An example of a Monte Carlo simulation is the one applied to approximate the value of pi. The simulation is based on generating random points within a unit square and see how many points fall within the circle enclosed by the unit square (marked in red). The higher the number of sampled points the closer the result is to the actual result. After selecting 30,000 random points, the estimate for pi is much closer to the actual value within the four decimal points of precision.
In this lab, we will learn how to generate random samples with various simulations and how to run a sensitivity analysis on the marketing use case covered so far.
Setup
Remember to always set your working directory to the source file location. Go to ‘Session’, scroll down to ‘Set Working Directory’, and click ‘To Source File Location’. Read carefully the below and follow the instructions to complete the tasks and answer any questions. Submit your work to RPubs as detailed in previous notes.
Note
For your assignment you may be using different data sets than what is included here. Always read carefully the instructions on Sakai. Tasks/questions to be completed/answered are highlighted in larger bolded fonts and numbered according to their particular placement in the task section.
PART A: SENSITIVITY ANALYSIS
In order to conduct the sensitivity analysis, we will need to download again the lpSolveAPI package unless you have it already installed in your R environment
# Require will load the package only if not installed
# Dependencies = TRUE makes sure that dependencies are install
if(!require("lpSolveAPI",quietly = TRUE))
install.packages("lpSolveAPI",dependencies = TRUE, repos = "https://cloud.r-project.org")We will revisit and solve again the marketing case discussed in class (also part of previous lab).
# We start with `0` constraint and `2` decision variables. The object name `lpmark` is discretionary.
lpmark = make.lp(0, 2)
# Define type of optimization as maximum and dump the screen output into a `dummy` variable
dummy = lp.control(lpmark, sense="max")
# Set the objective function coefficients
set.objfn(lpmark, c(275.691, 48.341))
#Add all constraints to the model
add.constraint(lpmark, c(1, 1), "<=", 350000)
add.constraint(lpmark, c(1, 0), ">=", 15000)
add.constraint(lpmark, c(0, 1), ">=", 75000)
add.constraint(lpmark, c(2, -1), "=", 0)
add.constraint(lpmark, c(1, 0), ">=", 0)
add.constraint(lpmark, c(0, 1), ">=", 0)
#Show the problem setting in tabular/matrix form. It's useful to see if our contraints have been properly set.
lpmark## Model name:
## C1 C2
## Maximize 275.691 48.341
## R1 1 1 <= 350000
## R2 1 0 >= 15000
## R3 0 1 >= 75000
## R4 2 -1 = 0
## R5 1 0 >= 0
## R6 0 1 >= 0
## Kind Std Std
## Type Real Real
## Upper Inf Inf
## Lower 0 0#Solve the linear programming problem
solve(lpmark)## [1] 0#The next two lines of codes will show the optimum results.
#Frist: Display the objective function optimum value i.e. the optimum sales value.
get.objective(lpmark)## [1] 43443517#Second: Display the decision variables optimum values i.e. the optimum values for radio and tv ads.
get.variables(lpmark) ## [1] 116666.7 233333.3For the sensitivity part we will add two new code sections to obtain the sensitivity results.
#Display sensitivity to the COEFFICIENTS of objective function.
get.sensitivity.obj(lpmark)## $objfrom
## [1] -96.6820 -137.8455
##
## $objtill
## [1] 1e+30 1e+30TASK 1: The results have two parts which are: the output labeled objfrom shows the lower limit of the coefficients while the output labeled objtill shows the upper limit. Explain in concise manner what the sensitivity results represent in reference to the marketing model.
ANSWER TASK 1:
The sensitivity results represent the lower limits of each of our coefficients, radio and tv. The lower limit for radio is -96.682 and the lower limit for tv is -137.846. They have no upper limits, which means that the values radio and tv can move anywhere in the scale from -96.682 to infinity or from -137.846 to infinity without changing the optimum solution. If they go down beyond the lower limit, the optimum solution will change. For every additional dollar spent on decision variable X1 more than the optimal solution, the optimal sales will decrease by about $96.68. For every additional dollar spent on decision variable X2 more than the optimal solution, the optimal sales will decrease by about $137.85.
#Display sensitivity to the CONSTRAINTS (or the right hand side values).
#There will be a total of m+n values where m is the number of contraints and n is the number of decision variables
get.sensitivity.rhs(lpmark) ## $duals
## [1] 124.12433 0.00000 0.00000 75.78333 0.00000 0.00000 0.00000
## [8] 0.00000
##
## $dualsfrom
## [1] 1.125e+05 -1.000e+30 -1.000e+30 -3.050e+05 -1.000e+30 -1.000e+30
## [7] -1.000e+30 -1.000e+30
##
## $dualstill
## [1] 1.00e+30 1.00e+30 1.00e+30 4.75e+05 1.00e+30 1.00e+30 1.00e+30 1.00e+30TASK 2: For this exercise we are only interested in the first part of the output which is labeled duals. Explain in coincise manner what the two non-zero sensitivity results represent. In your answer, distinguish between the binding and non-binding constraints, and include the explanation about the surplus/slack, and marginal values.
ANSWER TASK 2:
124.12433 and 75.78333 are the two non-zero sensitivity results which represent the binding conditions. Any change in those constraints will change the optimum value, which is highly sensitive to constraint changes. If budget is increased by $1, the optimum solution increases by $124.12433 (which is our marginal value). It is a binding constraint. To acquire a better understanding of the sensitivity results, and to confirm the integrity of the calculations, independent tests can be conducted. Binding constraints/resources are when the resource is limiting, whereby any increase will impact the optimal solution and there is no surplus in resource. Non-binding constraints/resources are when increasing the constraint has no impact on optimal solution and there is slack in resource.
TASK 3: Run the linear programing solver again starting from the begining, by defining a new model object lpmark1. All being equal, change the budget constraint by only $1 and solve. Specifially, all being equal, change the first constraint X1 + X2 <= 350000 by only $1 so that the new constraint will be X1 + X2 <= 350001. Note the optimum value for sales as given by the objective function.
# Define a new model object called lpmark1
lpmark1 = make.lp(0, 2)
# Repeat rest of commands with the one constraint change for budget. Solve and display the objective function optimum value
dummy = lp.control(lpmark1, sense="max")
set.objfn(lpmark1, c(275.691, 48.341))
add.constraint(lpmark1, c(1, 1), "<=", 350001)
add.constraint(lpmark1, c(1, 0), ">=", 15000)
add.constraint(lpmark1, c(0, 1), ">=", 75000)
add.constraint(lpmark1, c(2, -1), "=", 0)
add.constraint(lpmark1, c(1, 0), ">=", 0)
add.constraint(lpmark1, c(0, 1), ">=", 0)
lpmark1## Model name:
## C1 C2
## Maximize 275.691 48.341
## R1 1 1 <= 350001
## R2 1 0 >= 15000
## R3 0 1 >= 75000
## R4 2 -1 = 0
## R5 1 0 >= 0
## R6 0 1 >= 0
## Kind Std Std
## Type Real Real
## Upper Inf Inf
## Lower 0 0solve(lpmark1)## [1] 0get.objective(lpmark1)## [1] 43443641get.variables(lpmark1) ## [1] 116667 233334ANSWER TASK 4:
The new optimum value for sales is $4,344,361. The differential change in sales is $43,443,641 - $43,443,517, which is $124. That $124 is the same marginal value we thought we would get before using sensitivity analysis in TASK 2 located in $duals [1]. Changing the constraint changed our optimum value, which shows us that the optimum solution is sensitive to its binding conditions.
TASK 5: Running the linear programing solver again starting from the beginning, by defining a new model object lpmark2. All being equal, change the fourth constraint 2X1 - X2 = 0 by only $1 and solve. The new constraint will be 2X1 - X2 = 1. Note the optimum value for sales as given by the objective function.
# Define a new model object called lpmark2
lpmark2 = make.lp(0, 2)
# Repeat rest of commands with the above constraint changed. Solve and display the objective function optimum value
dummy = lp.control(lpmark2, sense="max")
set.objfn(lpmark2, c(275.691, 48.341))
add.constraint(lpmark2, c(1, 1), "<=", 350000)
add.constraint(lpmark2, c(1, 0), ">=", 15000)
add.constraint(lpmark2, c(0, 1), ">=", 75000)
add.constraint(lpmark2, c(2, -1), "=", 1)
add.constraint(lpmark2, c(1, 0), ">=", 0)
add.constraint(lpmark2, c(0, 1), ">=", 0)
lpmark2## Model name:
## C1 C2
## Maximize 275.691 48.341
## R1 1 1 <= 350000
## R2 1 0 >= 15000
## R3 0 1 >= 75000
## R4 2 -1 = 1
## R5 1 0 >= 0
## R6 0 1 >= 0
## Kind Std Std
## Type Real Real
## Upper Inf Inf
## Lower 0 0solve(lpmark2)## [1] 0get.objective(lpmark2)## [1] 43443592get.variables(lpmark2)## [1] 116667 233333ANSWER TASK 6:
The new optimum value for sales is $4,344,361. The differential change in sales is $43,443,592 - $43,443,517, which is $75. That $75 is the same marginal value we thought we would get before using sensitivity analysis in TASK 2 located in $duals [1]. Changing the constraint changed our optimum value, which shows us that the optimum solution is not sensitive to its conditions and that the conditions are non-binding. When Constraint R4 is increased by 1, the optimum sales value increases by $75.
PART B: MONTE CARLO SIMULATION
For this task we will be running a Monte Carlo simulation to calculate the probability that the daily return from S&P will be > 5%. We will assume that the historical S&P daily return follows a normal distribution with an average daily return of 0.03 (%) and a standard deviation of 0.97 (%).
To begin we will generate 100 random samples from the normal distribution. For the generated samples we will calculate the mean, standard deviation, and probability of occurrence where the simulation result is greater than 5%.
To generate random samples from a normal distribution we will use the rnorm() function in R. In the example below we set the number of runs (or samples) to 100.
# number of simulations/samples
runs = 100
# random number generator per defined normal distribution with given mean and standard deviation
sims = rnorm(runs,mean=0.03,sd=0.97)# Mean calculated from the random distribution of samples
average = mean(sims)
average## [1] 0.06739668# STD calculated from the random distribution of samples
std = sd(sims)
std## [1] 1.062748# probability of occurrence on any given day based on samples will be equal to count (or sum) where sample result is greater than 5% divided by total number of samples.
prob = sum(sims >=0.05)/runs
prob## [1] 0.45TASK 7: Repeat the above calculations for the case where the number of simulations/samples is equal to 1000. record the mean, standard deviation, and probability. Name all the required variables as runs1, sims1, average1, std1, and prob1.
# Repeat calculations here
runs1 = 1000
sims1 = rnorm(runs,mean=0.03,sd=0.97)
average1 = mean(sims1)
average1## [1] -0.003137732std1 = sd(sims1)
std1## [1] 0.8633306prob1 = sum(sims1 >=0.05)/runs
prob1## [1] 0.47TASK 8: Repeat the above calculations for the case where the number of simulations/samples is equal to 10000. record the mean, standard deviation, and probability. Name all the required variables as runs2, sims2, average2, std2, and prob2.
# Repeat calculations here
runs2 = 10000
sims2 = rnorm(runs,mean=0.03,sd=0.97)
average2 = mean(sims2)
average2## [1] 0.2113135std2 = sd(sims2)
std2## [1] 0.8960543prob2 = sum(sims2 >=0.05)/runs
prob2## [1] 0.61TASK 9: List in a tabular form the values for mean, standard deviation, and probability for all three cases: 100, 1000, and 10000 simulations.
Runs = c("100","1000","10000")
Average = c(average,average1,average2)
STD = c(std,std1,std2)
Probability = c(prob,prob1,prob2)
myresults = data.frame(Runs,Average,STD,Probability)
myresults## Runs Average STD Probability
## 1 100 0.067396683 1.0627484 0.45
## 2 1000 -0.003137732 0.8633306 0.47
## 3 10000 0.211313463 0.8960543 0.61TASK 10: Describe how the values change/behave as the number of simulations is increased. What is your best bet on the probability of occurrence greater than 5% and why? How is this similar to the image use case to calculate pi that was presented in the introductory paragraph?
ANSWER TASK 10:
The probability of occurrence is always between [0,1]. The more scenarios we pick, the more accurate our computed average from the sample will be to the real population average and average standard deviation. The probability of occurrence being greater than 5% is fairly likely. The pi image shows a maximum of 1.0 on either axis, and as we increase the number of samples, we are likely to get closer to being very certainly having a probability of occurrence greater than 5%. The image used to calculate pi was stated, “The higher the number of sampled points, the closer the result is to the actual result. After selecting 30,000 random points, the estimate for pi is much closer to the actual value within the four decimal points of precision.” This is true for our case as well. After selecting 10,000 random points, our estimate would likely be more accurate.
The last 2C) exercise is optional for those interested in further enhancing their subject matter learning, and refining their skills in R. Your work will be assessed but you will not be graded for this exercise. You can follow the instructions presented in the video Excel equivalent example at [https://www.youtube.com/watch?v=wKdmEXCvo9s]