(homework) Decision Science #3 - NCCU (Howard Chuang)
JY
11/23/2020
Q1
Let’s revisit the case of Conley Fisheries. Define DGlou~triangular(min=4000, max=8000, mode=7000) and DRock~triangular (min=4800, max=7200, mode=6300) as stochastic demand for cold fish. The prices are random PriceGlou~Normal(μ=3.5/lb, σ=0.35/lb) and PriceRock~Normal(μ=3.65/lb, σ=0.25/lb). Also, corr(PriceGlou, PriceRock) is not zero. In addition, the Conley Fisheries has two boats (each with a full load of 3800 lbs) driven by captain Rick and captain Morty. The daily fraction captured by Risk is a uniform(0.6, 0.9) random variable whereas the daily fraction captured by Morty is a triangular(min=0.5, max=1, mode=0.75) random variable. The day-to-day operating costs for each boat is about $7,200.
Being the business analyst of Conley Fisheries, you lay out four possible selling strategies: - a) Rick goes to Glou and Morty goes to Rock; - b) Rick goes to Rock and Morty goes to Glou; - c) Rick and Morty go to Glou; - d) Rick and Morty go to Rock.
To assess which strategy is more profitable/less risky, write a simulation program that allows one to compute expected profit and CVAR(5%) for the four selling strategies above. After that, please finish the following tasks.
import packages
# install.packages('EnvStats')
# install.packages('MASS')
# install.packages('RcmdrMisc')
library(EnvStats)
library(MASS)
library(RcmdrMisc)set variables
### S: The number of simulation runs
S = 10000
output_option = 'mu'
strategy = 'a'
corrPGxPR = seq(-0.8, 0.8, by=0.2)define functions
mu.profit <- function(S, corrPGxPR, strategy, output_option) {
### market info
demand.Glou=round(rtri(S,4000,8000,7000),0)
mu.price.Glou = 3.5; sigma.price.Glou = 0.35
#sim.price.Glou = rnorm(S, mu.price.Glou, sigma.price.Glou)
demand.Rock=round(rtri(S,4800,7200,6300),0)
mu.price.Rock = 3.65; sigma.price.Rock = 0.25
#sim.price.Rock = rnorm(S, mu.price.Rock, sigma.price.Rock)
corrPGxPR = corrPGxPR
covPGxPR = sigma.price.Glou*sigma.price.Rock*corrPGxPR
covMatrix = matrix(c(sigma.price.Glou^2, covPGxPR, covPGxPR, sigma.price.Rock^2), nrow=2)
sim.price = mvrnorm(S, c(mu.price.Glou, mu.price.Rock), covMatrix)
sim.price.Glou = sim.price[,1]; sim.price.Rock = sim.price[,2]
#x11(width=12,height=5)
#par(mfrow=c(1,2))
#hist(demand.Glou)
#hist(demand.Rock)
### fish catching capacity
fullload.Rick = 3800
frac.Rick=runif(S, 0.6, 0.9)
fish.Rick=round(fullload.Rick*frac.Rick, 0)
fullload.Morty = 3800
frac.Morty=rtri(S, 0.5, 1, 0.75)
fish.Morty=round(fullload.Morty*frac.Morty, 0)
operationcost = 7200
if(strategy == 'a') {
### Situation (a)
Rick.Glou=c() # Rick.Glou is the stochastic revenue when selling all the fish at Glouport
Morty.Rock=c() # Morty.Rock is the stochastic revenue when selling all the fish at Rockport
for(s in 1:S){
Rick.Glou[s] = sim.price.Glou[s]*min(fish.Rick[s], demand.Glou[s]) - operationcost
Morty.Rock[s] = sim.price.Rock[s]*min(fish.Morty[s], demand.Rock[s]) - operationcost
}
summary(Rick.Glou); summary(Morty.Rock)
#x11(width=12,height=5)
#par(mfrow=c(1,2))
#hist(Rick.Glou, breaks=200)
#hist(Morty.Rock, breaks=200)
total.profit = Rick.Glou+Morty.Rock
} else if(strategy == 'b') {
### Situation (b)
Rick.Rock=c() # Rick.Rock is the stochastic revenue when selling all the fish at Rockport
Morty.Glou=c() # Morty.Glou is the stochastic revenue when selling all the fish at Glouport
for(s in 1:S){
Rick.Rock[s] = sim.price.Rock[s]*min(fish.Rick[s], demand.Rock[s]) - operationcost
Morty.Glou[s] = sim.price.Glou[s]*min(fish.Morty[s], demand.Glou[s]) - operationcost
}
summary(Rick.Rock); summary(Morty.Glou)
total.profit = Rick.Rock+Morty.Glou
} else if(strategy == 'c') {
### Situation (c)
Rick.Glou=c() # Rick.Glou is the stochastic revenue when selling all the fish at Glouport
Morty.Glou=c() # Morty.Glou is the stochastic revenue when selling all the fish at Glouport
for(s in 1:S){
Rick.Glou[s] = sim.price.Glou[s]*min(fish.Rick[s], demand.Glou[s]) - operationcost
Morty.Glou[s] = sim.price.Glou[s]*min(fish.Morty[s], demand.Glou[s]-min(fish.Rick[s], demand.Glou[s])) - operationcost
}
summary(Rick.Glou); summary(Morty.Glou)
total.profit = Rick.Glou+Morty.Glou
} else {
### Situation (d)
Rick.Rock=c() # Rick.Rock is the stochastic revenue when selling all the fish at Rockport
Morty.Rock=c() # Morty.Rock is the stochastic revenue when selling all the fish at Rockport
for(s in 1:S){
Rick.Rock[s] = sim.price.Rock[s]*min(fish.Rick[s], demand.Rock[s]) - operationcost
Morty.Rock[s] = sim.price.Rock[s]*min(fish.Morty[s], demand.Rock[s]-min(fish.Rick[s], demand.Rock[s])) - operationcost
}
summary(Rick.Rock); summary(Morty.Rock)
total.profit = Rick.Rock+Morty.Rock
}
CVARq5 = mean(total.profit[which(total.profit <= quantile(total.profit, 0.05))])
if (output_option == 'mu') {
return(c(sum(total.profit)/S, cor(sim.price.Glou, sim.price.Rock)))
} else {
return(c(CVARq5, cor(sim.price.Glou, sim.price.Rock)))
}
}Situation (a)
a.muprofit = c(); a.cvarq5 = c()
for(s in 1:length(corrPGxPR)) {
a.muprofit[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='a', output_option='mu')
a.cvarq5[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='a', output_option='cvarq5')
}
cat("average of profit:", mean(a.muprofit), "\n")## average of profit: 5976.812cat("CVAR(q=5%) :", mean(a.cvarq5), "\n")## CVAR(q=5%) : 1686.107Situation (b)
b.muprofit = c(); b.cvarq5 = c()
for(s in 1:length(corrPGxPR)) {
b.muprofit[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='b', output_option='mu')
b.cvarq5[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='b', output_option='cvarq5')
}
cat("average of profit:", mean(b.muprofit), "\n")## average of profit: 5980.913cat("CVAR(q=5%) :", mean(b.cvarq5), "\n")## CVAR(q=5%) : 1734.586Situation (c)
c.muprofit = c(); c.cvarq5 = c()
for(s in 1:length(corrPGxPR)) {
c.muprofit[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='c', output_option='mu')
c.cvarq5[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='c', output_option='cvarq5')
}
cat("average of profit:", mean(c.muprofit), "\n")## average of profit: 4950.089cat("CVAR(q=5%) :", mean(c.cvarq5), "\n")## CVAR(q=5%) : -238.5945Situation (d)
d.muprofit = c(); d.cvarq5 = c()
for(s in 1:length(corrPGxPR)) {
d.muprofit[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='d', output_option='mu')
d.cvarq5[s] = mu.profit(S, corrPGxPR=corrPGxPR[s], strategy='d', output_option='cvarq5')
}
cat("average of profit:", mean(d.muprofit), "\n")## average of profit: 5923.478cat("CVAR(q=5%) :", mean(d.cvarq5), "\n")## CVAR(q=5%) : 1787.158Q1-1
- (20%) Show a figure where the y-axis is expected profit and x-axis is corr(PriceGlou, PriceRock). For each of the following values in the x-axis: -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, and 0.8, simulate 10,000 runs respectively. Compute and plot the expected profits of the four strategies and. The figure should show four lines (in different colors and types).
lineplot(corrPGxPR, a.muprofit, b.muprofit, c.muprofit, d.muprofit)
Q1-2
- (20%) Show a figure where the y-axis is CVAR(5%) and x-axis is corr(PriceGlou, PriceRock). For each of the following values in the x-axis: -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, and 0.8, simulate 10,000 runs respectively. Compute and plot CVAR(5%) values of the four strategies. The figure should exhibit four lines (in different colors and types).
lineplot(corrPGxPR, a.cvarq5, b.cvarq5, c.cvarq5, d.cvarq5)
Q1-3
- (10%) Provide a succinct and logical discussion about findings in the two figures above. Also, explain what kind of selling strategies you would recommend the CEO to adopt.
For more risky strategy, observe the expected value of profit, except for situation (c), other strategy can gain higher average profit;
For more conservative strategy, observe the CVAR(q=5%),
when the corr(PriceGlou, PriceRock) getting higher, choose situation (d) will get higher CVAR(q=5%), when the corr(PriceGlou, PriceRock) getting lower, choose either situation (a) or situation (b) will get higher CVAR(q=5%).
—
Q2
NCCU Griffins (雄鷹) will play an important Bowl game next month, and we have an official license to sell official shirts. The shirts cost $10 each to produce, and we will sell them before the Bowl game at a price of $25. At this price, the anticipated demand before the game should be a Normal(9000, 2000) random variable (RV). After the Bowl game, this depend on whether NCCU Griffins win or not. If they win, we will continue to sell the shirt for $25, and demand in the month afterwards will be a Normal(6000, 2000) RV. If they lose, we will cut the price to $12.5 and the demand should be a gamma RV with mean 2000 and standard deviation of 1000. The probability of our local heroes winning is 0.4. Unsold shirts will be given away for free.
set variables
#cost, price
shirt.cost = 10; shirt.price = 25; race.lose.shirt.price = shirt.price/2
#demend
shirt.demand.before = round(rnorm(n=S, mean=9000, sd=2000), 0); shirt.demand.before[shirt.demand.before<0] <- 0
race.win.shirt.demand = round(rnorm(n=S, mean=6000, sd=2000), 0); race.win.shirt.demand[race.win.shirt.demand<0] <- 0
race.lose.shirt.demand = round(rnorm(n=S, mean=2000, sd=1000), 0); race.lose.shirt.demand[race.lose.shirt.demand<0] <- 0
#chance to win
perfectinfo = TRUE
race.win.perfect = 0.4
race.win = runif(n=1, min=0, max=1)
winorlose.perfect = sample(x=c(1,0), size=S, prob=c(race.win.perfect, 1-race.win.perfect), replace=TRUE)
winorlose = sample(x=c(1,0), size=S, prob=c(race.win, 1-race.win), replace=TRUE)
#production
shirt.produce = round(rtri(
20,
max=round(max(shirt.demand.before+race.win.shirt.demand), 0),
min=round(min(shirt.demand.before+race.lose.shirt.demand), 0),
mode=round(mean(shirt.demand.before+(0.4*race.win.shirt.demand+0.6*race.lose.shirt.demand)), 0)
), 0)
shirt.produce = seq(round(min(shirt.demand.before+race.lose.shirt.demand), 0),
round(max(shirt.demand.before+race.win.shirt.demand), 0),
by=1000)define functions
### demand function
demand.perfect <- function(S, winorlose.perfect, shirt.demand.before, race.win.shirt.demand, race.lose.shirt.demand){
demand = c()
for (s in 1:S) {
if (winorlose.perfect[s] == 1) {
demand[s] = shirt.demand.before[s] + race.win.shirt.demand[s]
} else {
demand[s] = shirt.demand.before[s] + race.lose.shirt.demand[s]
}
}
return(demand)
}
### profit function
profit <- function(shirt.cost, shirt.produce,
shirt.price, race.lose.shirt.price,
winorlose, perfectinfo,
shirt.demand.before, race.win.shirt.demand, race.lose.shirt.demand){
sim.profit.before = c()
sim.profit.after = c()
for(s in 1:S) {
sim.profit.before[s] = (shirt.price-shirt.cost)*min(shirt.produce, shirt.demand.before[s])
shirt.left = shirt.produce - min(shirt.produce, shirt.demand.before[s])
if (perfectinfo == TRUE) { winorlose = winorlose.perfect }
if (winorlose[s] == 1) {
sim.profit.after[s] = (shirt.price-shirt.cost)*min(shirt.left, race.win.shirt.demand[s]) - (shirt.cost*max(0, shirt.left-race.win.shirt.demand[s]))
} else {
sim.profit.after[s] = (race.lose.shirt.price-shirt.cost)*min(shirt.left, race.lose.shirt.demand[s]) - (shirt.cost*max(0, shirt.left-race.lose.shirt.demand[s]))
}
}
return(sim.profit.before+sim.profit.after)
}start simulation
sim.profit = matrix(0, nrow=S, ncol=length(shirt.produce))
avg.profit = c()
sd.profit = c()
CVARq10=c()
for (i in 1:length(shirt.produce)) {
sim.profit[, i] = profit(shirt.cost, shirt.produce[i],
shirt.price, race.lose.shirt.price,
winorlose.perfect, perfectinfo,
shirt.demand.before, race.win.shirt.demand, race.lose.shirt.demand)
avg.profit[i] = mean(sim.profit[, i])
sd.profit[i] = sd(sim.profit[, i])
CVARq10[i] = mean(sim.profit[, i][which(sim.profit[, i] <= quantile(sim.profit[, i], 0.1))])
cat("production quantity:", shirt.produce[i], "\n")
}## production quantity: 2322
## production quantity: 3322
## production quantity: 4322
## production quantity: 5322
## production quantity: 6322
## production quantity: 7322
## production quantity: 8322
## production quantity: 9322
## production quantity: 10322
## production quantity: 11322
## production quantity: 12322
## production quantity: 13322
## production quantity: 14322
## production quantity: 15322
## production quantity: 16322
## production quantity: 17322
## production quantity: 18322
## production quantity: 19322
## production quantity: 20322
## production quantity: 21322
## production quantity: 22322
## production quantity: 23322
## production quantity: 24322
## production quantity: 25322lineplot(shirt.produce, avg.profit, sd.profit, CVARq10)
Q2-1
- (20%) What is the production quantity that would maximize our expected profit? Please define a reasonable range of possible quantities for search.
cat("the optimal production quantity for the maximize expected profit:", shirt.produce[which.max(avg.profit)], "\n")## the optimal production quantity for the maximize expected profit: 12322Q2-2
- (20%) How would the optimal production quantity change if we instead want to maximize CVAR(q=10%)?
cat("the optimal production quantity for the maximize CVAR(q=10%):", shirt.produce[which.max(CVARq10)], "\n")## the optimal production quantity for the maximize CVAR(q=10%): 6322Q2-3
- (10%) Assuming we want to maximize expected profit in (a), what would be the expected value of perfect information about whether our local heroes will win the Bowl or not?
### Assess the expected value perfect information
#### 先求出至少滿足服務水準，所需生產的數量
cu.before = shirt.price-shirt.cost
co.before = shirt.cost-0
frac.before = cu.before/(cu.before+co.before)
cu.lose = race.lose.shirt.price-shirt.cost
co.lose = shirt.cost-0
frac.lose = cu.lose/(cu.lose+co.lose)
stdlevel.produce.win = round(quantile(shirt.demand.before,frac.before,names=FALSE)) + round(quantile(race.win.shirt.demand,frac.before,names=FALSE), 0)
stdlevel.produce.lose = round(quantile(shirt.demand.before,frac.before,names=FALSE)) + round(quantile(race.lose.shirt.demand,frac.lose,names=FALSE), 0)
stdlevel.profit.win = profit(shirt.cost, stdlevel.produce.win,
shirt.price, race.lose.shirt.price,
winorlose, perfectinfo,
shirt.demand.before, race.win.shirt.demand, race.lose.shirt.demand)
mean(stdlevel.profit.win)## [1] 133150.4stdlevel.profit.lose = profit(shirt.cost, stdlevel.produce.lose,
shirt.price, race.lose.shirt.price,
winorlose, perfectinfo,
shirt.demand.before, race.win.shirt.demand, race.lose.shirt.demand)
mean(stdlevel.profit.lose)## [1] 139661.9#### 模擬在 win, lose 各自情況下，會選擇的策略和所得收益
profit.perfectinfo=c()
demand.perfect = demand.perfect(S, winorlose.perfect, shirt.demand.before, race.win.shirt.demand, race.lose.shirt.demand)
for(s in 1:S) {
if (winorlose.perfect[s] == 1) {
shirt.produce = stdlevel.produce.win
profit.before = (shirt.price-shirt.cost)*min(shirt.produce, shirt.demand.before[s])
shirt.left = shirt.produce - min(shirt.produce, shirt.demand.before[s])
profit.after = (shirt.price-shirt.cost)*min(shirt.left, race.win.shirt.demand[s]) - (shirt.cost*max(0, shirt.left-race.win.shirt.demand[s]))
} else {
shirt.produce = stdlevel.produce.lose
profit.before = (shirt.price-shirt.cost)*min(shirt.produce, shirt.demand.before[s])
shirt.left = shirt.produce - min(shirt.produce, shirt.demand.before[s])
profit.after = (race.lose.shirt.price-shirt.cost)*min(shirt.left, race.lose.shirt.demand[s]) - (shirt.cost*max(0, shirt.left-race.lose.shirt.demand[s]))
}
profit.perfectinfo[s] = profit.before + profit.after
}
cat("expected value of perfect information:", mean(profit.perfectinfo), "\n")## expected value of perfect information: 155444.7cat("maximum average profit:", max(avg.profit), "\n")## maximum average profit: 144400.8cat("how much we benefit from perfect information:", mean(profit.perfectinfo)-max(avg.profit), "\n")## how much we benefit from perfect information: 11043.88