Homework III (Group)

Due date: 23:59 on November 27 (Friday), 2020

library(EnvStats)

## 
## Attaching package: 'EnvStats'

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

## The following object is masked from 'package:base':
## 
##     print.default

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:EnvStats':
## 
##     boxcox

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Q1 (50%)

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 Rick 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.

S = 10000
fullload = 3800
ope.costs = 7200
# simulate demand for 2 ports
sim.D.Glou = round(rtri(S, min=4000, max=8000, mode=7000),0)
sim.D.Rock = round(rtri(S, min=4800, max=7200, mode=6300),0)

# simulate daily captured fish for 2 boats
frac.R = runif(S, 0.6, 0.9)
frac.M = rtri(S, min=0.5, max=1, mode=0.75)
sim.fish.R = round(fullload*frac.R,0)
sim.fish.M = round(fullload*frac.M,0)
# PriceGlou~Normal(μ=3.5/lb, σ=0.35/lb) and PriceRock~Normal(μ=3.65/lb, σ=0.25/lb)
mu.PRG = 3.5
mu.PRR = 3.65
sigma.PRG = 0.35
sigma.PRR = 0.25

corr.GR = seq(-0.8,0.8,0.2)

muRev = matrix(nrow=length(corr.GR),ncol=4)
CVaRq5 = matrix(nrow=length(corr.GR),ncol=4)
for (i in 1:length(corr.GR)) {
    # simulate different correlation and daily price
    covGR = sigma.PRG*sigma.PRR*corr.GR[i]
    cov.matrix = matrix(c(sigma.PRG^2, covGR, covGR, sigma.PRR^2),nrow=2)
    sim.PR = mvrnorm(S, c(mu.PRG,mu.PRR), cov.matrix)
    earnings = matrix(nrow=S,ncol=4)
    # simulate 4 strategies
    for (s in 1:S) {
        # Rick goes to Glou and Morty goes to Rock
        earnings[s,1]=sim.PR[s,1]*min(sim.D.Glou[s], sim.fish.R[s]) + sim.PR[s,2]*min(sim.D.Rock[s],sim.fish.M[s])-2*ope.costs
        # Rick goes to Rock and Morty goes to Glou
        earnings[s,2]=sim.PR[s,2]*min(sim.D.Rock[s], sim.fish.R[s]) + sim.PR[s,1]*min(sim.D.Glou[s],sim.fish.M[s])-2*ope.costs
        # Rick and Morty go to Glou
        earnings[s,3]=sim.PR[s,1]*min(sim.D.Glou[s], (sim.fish.M[s]+sim.fish.R[s]))-2*ope.costs
        # Rick and Morty go to Rock
        earnings[s,4]=sim.PR[s,2]*min(sim.D.Rock[s], (sim.fish.M[s]+sim.fish.R[s]))-2*ope.costs
    }
    # for each strategy i: correlation
    for (j in 1:4) {
        muRev[i,j] = mean(earnings[,j])
        valq5 = quantile(earnings[,j],0.05)
        CVaRq5[i,j] = mean(earnings[which(earnings[,j]<=valq5),j])
    }
    cat("corr=",corr.GR[i], "done\n")
}

## corr= -0.8 done
## corr= -0.6 done
## corr= -0.4 done
## corr= -0.2 done
## corr= 0 done
## corr= 0.2 done
## corr= 0.4 done
## corr= 0.6 done
## corr= 0.8 done

(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).

plot(corr.GR, muRev[,1], "l",col="red", ylim=c(4500, 6500),xaxt="n")
lines(corr.GR, muRev[,2],col="blue")
lines(corr.GR, muRev[,3],col="green")
lines(corr.GR, muRev[,4],col="black")
axis(1, at = corr.GR, las=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).

plot(corr.GR, CVaRq5[,1], "l",col="red", ylim=c(-4000, 4000),xaxt="n")
lines(corr.GR, CVaRq5[,2],col="blue")
lines(corr.GR, CVaRq5[,3],col="green")
lines(corr.GR, CVaRq5[,4],col="black")
axis(1, at = corr.GR, las=2)

(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.

A: 兩條船要各自去不同的港口才會賺錢，不要去同一個港口，至於誰去哪個港口都沒差

Q2 (50%)

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.(salvage value=0)

\[ Maximize \; (price-cost)*min(D,Q)-(c-sv)*max(Q-D,0) \]

S = 10000
scale = 1000^2/2000
shape = 2000/scale
sim.D.before = round(rnorm(S, 9000,2000),0)
# prob(win)
p.win = 0.4

# 1 denotes win; 0 denotes lose
sim.bowlgame = sample(c(1,0), S, prob=c(p.win, 1-p.win), replace=T)
sim.D.after = matrix(0, nrow=S, ncol=2)
colnames(sim.D.after) = c("isWin","demand")
for (i in 1:S) {
    # win
    if (sim.bowlgame[i] == 1) {
        sim.D.after[i,1] = 1
        sim.D.after[i,2] = rnorm(1, 6000,2000)
    }
    else {
        sim.D.after[i, 1] = 0
        sim.D.after[i, 2]=rgamma(1, shape=shape, scale=scale)
    }
    sim.D.after[i,2]=round(sim.D.after[i,2], 0)
}

unitcost = 10
unitprice = 25
salvage = 0
F.profit = function(Q, D1, D2) {
    profit.val=c()
    for(i in 1:length(D1)) {
        soldQ = min(Q, D1[i])
        revenue1=(unitprice-unitcost)*soldQ
        remainQ=Q-soldQ
        if (D2[i, "isWin"]==1) {
            revenue2=(unitprice-unitcost)*min(remainQ, D2[i, "demand"])+(salvage-unitcost)*max(remainQ-D2[i, "demand"], 0)
        }
        else {
            revenue2=(unitprice/2-unitcost)*min(remainQ, D2[i, "demand"])+(salvage-unitcost)*max(remainQ-D2[i, "demand"], 0)
        }
        profit.val[i]= revenue1 + revenue2
    }
    return(profit.val)
}

Q.val=seq(3000,100000,1000)
sim.profit=matrix(0,nrow=S,ncol=length(Q.val))
colnames(sim.profit)=Q.val
avg.profit=c()
sd.profit=c()

start.time <- Sys.time()

for(i in 1:length(Q.val)){
  sim.profit[,i]=F.profit(Q.val[i], sim.D.before, sim.D.after)
  avg.profit[i]=mean(sim.profit[,i])
  sd.profit[i]=sd(sim.profit[,i])
  cat("production quantity:", Q.val[i], "\n")
}

## production quantity: 3000 
## production quantity: 4000 
## production quantity: 5000 
## production quantity: 6000 
## production quantity: 7000 
## production quantity: 8000 
## production quantity: 9000 
## production quantity: 10000 
## production quantity: 11000 
## production quantity: 12000 
## production quantity: 13000 
## production quantity: 14000 
## production quantity: 15000 
## production quantity: 16000 
## production quantity: 17000 
## production quantity: 18000 
## production quantity: 19000 
## production quantity: 20000 
## production quantity: 21000 
## production quantity: 22000 
## production quantity: 23000 
## production quantity: 24000 
## production quantity: 25000 
## production quantity: 26000 
## production quantity: 27000 
## production quantity: 28000 
## production quantity: 29000 
## production quantity: 30000 
## production quantity: 31000 
## production quantity: 32000 
## production quantity: 33000 
## production quantity: 34000 
## production quantity: 35000 
## production quantity: 36000 
## production quantity: 37000 
## production quantity: 38000 
## production quantity: 39000 
## production quantity: 40000 
## production quantity: 41000 
## production quantity: 42000 
## production quantity: 43000 
## production quantity: 44000 
## production quantity: 45000 
## production quantity: 46000 
## production quantity: 47000 
## production quantity: 48000 
## production quantity: 49000 
## production quantity: 50000 
## production quantity: 51000 
## production quantity: 52000 
## production quantity: 53000 
## production quantity: 54000 
## production quantity: 55000 
## production quantity: 56000 
## production quantity: 57000 
## production quantity: 58000 
## production quantity: 59000 
## production quantity: 60000 
## production quantity: 61000 
## production quantity: 62000 
## production quantity: 63000 
## production quantity: 64000 
## production quantity: 65000 
## production quantity: 66000 
## production quantity: 67000 
## production quantity: 68000 
## production quantity: 69000 
## production quantity: 70000 
## production quantity: 71000 
## production quantity: 72000 
## production quantity: 73000 
## production quantity: 74000 
## production quantity: 75000 
## production quantity: 76000 
## production quantity: 77000 
## production quantity: 78000 
## production quantity: 79000 
## production quantity: 80000 
## production quantity: 81000 
## production quantity: 82000 
## production quantity: 83000 
## production quantity: 84000 
## production quantity: 85000 
## production quantity: 86000 
## production quantity: 87000 
## production quantity: 88000 
## production quantity: 89000 
## production quantity: 90000 
## production quantity: 91000 
## production quantity: 92000 
## production quantity: 93000 
## production quantity: 94000 
## production quantity: 95000 
## production quantity: 96000 
## production quantity: 97000 
## production quantity: 98000 
## production quantity: 99000 
## production quantity: 1e+05

end.time <- Sys.time()
cat(end.time-start.time)

## 4.998447

par(mfrow=c(2,1))
plot(Q.val, avg.profit,type='l',xlab="production quantity",lwd=3)
# best point
points(Q.val[which.max(avg.profit)], max(avg.profit), col="red",pch=19,cex=1.2)
# zoom in
plot(Q.val[which(avg.profit>0)], avg.profit[which(avg.profit>0)], type='l',xlab="production quantity",lwd=3)

Q.val[which(avg.profit>0)]

##  [1]  3000  4000  5000  6000  7000  8000  9000 10000 11000 12000 13000 14000
## [13] 15000 16000 17000 18000 19000 20000 21000 22000 23000 24000 25000 26000
## [25] 27000 28000 29000 30000

可以看出主要賺錢的範圍為生產3000~29000，期望利潤最高的為生產13000，我們可以繼續往下搜尋更精密的數量

Q.val=seq(3000,29000,100)
sim.profit=matrix(0,nrow=S,ncol=length(Q.val))
colnames(sim.profit)=Q.val
avg.profit=c()
sd.profit=c()
CVaRq10=c()
start.time <- Sys.time()

for(i in 1:length(Q.val)){
  sim.profit[,i]=F.profit(Q.val[i], sim.D.before, sim.D.after)
  avg.profit[i]=mean(sim.profit[,i])
  sd.profit[i]=sd(sim.profit[,i])
  lowestq10 = sim.profit[which(sim.profit[,i]<=quantile(sim.profit[,i], 0.1)),i]
  CVaRq10[i] = mean(lowestq10)
  if (Q.val[i]%%1000==0) {
    cat("production quantity:", Q.val[i], "\n")
  }
}

## production quantity: 3000 
## production quantity: 4000 
## production quantity: 5000 
## production quantity: 6000 
## production quantity: 7000 
## production quantity: 8000 
## production quantity: 9000 
## production quantity: 10000 
## production quantity: 11000 
## production quantity: 12000 
## production quantity: 13000 
## production quantity: 14000 
## production quantity: 15000 
## production quantity: 16000 
## production quantity: 17000 
## production quantity: 18000 
## production quantity: 19000 
## production quantity: 20000 
## production quantity: 21000 
## production quantity: 22000 
## production quantity: 23000 
## production quantity: 24000 
## production quantity: 25000 
## production quantity: 26000 
## production quantity: 27000 
## production quantity: 28000 
## production quantity: 29000

end.time <- Sys.time()
cat(end.time-start.time)

## 13.17084

plot(Q.val, avg.profit,type='l',xlab="production quantity",lwd=3)
# best point
points(Q.val[which.max(avg.profit)], max(avg.profit), col="red",pch=19,cex=2)

(20%) What is the production quantity that would maximize our expected profit? Please define a reasonable range of possible quantities for search.

cat("Produce quantity of",Q.val[which.max(avg.profit)], " E(x):", max(avg.profit))

## Produce quantity of 12600  E(x): 144597.8

Produce 12,600 will be the maximum of expected profit.

(20%) How would the optimal production quantity change if we instead want to maximize CVAR(q=10%)?

cat("Produce at risk 10% quantity of",Q.val[which.max(CVaRq10)], " E(x):", max(CVaRq10))

## Produce at risk 10% quantity of 7000  E(x): 103570.4

plot(Q.val, CVaRq10,type='l',xlab="production quantity",lwd=3)
# best point
points(Q.val[which.max(CVaRq10)], max(CVaRq10), col="red",pch=19,cex=2)

(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?

事前知道雄鷹可能會贏或是輸

cu=unitprice-unitcost
co=unitcost-salvage
frac=cu/(cu+co)
cat(frac)

## 0.6

cu2=unitprice/2-unitcost
frac2=cu2/(cu2+co)
cat(frac2)

## 0.2

x.before = quantile(sim.D.before, frac,names=F) %>% round(0)
x.after.win = quantile(sim.D.after[sim.D.after[,1]==1,"demand"], frac,names=F) %>% round(0)
x.after.lose = quantile(sim.D.after[sim.D.after[,1]==0,"demand"], frac2,names=F) %>% round(0)

profit.perfectinfo=c()

for(i in 1:length(sim.D.before)) {
    soldQ = min(x.before, sim.D.before[i])
    revenue1=(unitprice-unitcost)*soldQ
    if (sim.D.after[i, "isWin"]==1) {
        revenue2=(unitprice-unitcost)*min(x.after.win, sim.D.after[i, "demand"])+
          (salvage-unitcost)*max(x.after.win-sim.D.after[i, "demand"], 0)
    }
    else {
        revenue2=(unitprice/2-unitcost)*min(x.after.lose, sim.D.after[i, "demand"])+
          (salvage-unitcost)*max(x.after.lose-sim.D.after[i, "demand"], 0)
    }
    profit.perfectinfo[i]= revenue1 + revenue2
}

mean(profit.perfectinfo)-max(avg.profit)

## [1] 11278.74

1091 NCCU Decision Science HW3

Homework III (Group)

Q1 (50%)

Q2 (50%)