# 決策科學HW3
## Q1

%pip install scipy --upgrade
import math
import matplotlib.pyplot as plt
import scipy.stats as stats
import scipy.optimize as opt
import numpy as np
import pandas as pd
import statsmodels.api as sm
%pip install pymoo
import pymoo as pymoo

#define a demand function
a=200
b=35

def Demand(Price,a=200,b=35):
  Dval=a-b*Price+np.random.normal(0,20)
  Dval=np.round(Dval)
  return max(Dval,0)

S=5000
salvage_value=0.5
shortage_cost=1
purchase_cost=1
def Perf(x):
  np.random.seed(5566)
  P=x[0]
  Q=x[1]
  profit=np.empty(S)
  for s in range(S):
    demand_s=Demand(P)
    profit[s]=P*min(Q,demand_s) - purchase_cost*Q + salvage_value*max(Q-demand_s,0) - shortage_cost*max(demand_s-Q,0)
  avgprofit=np.mean(profit)
  return -1*avgprofit

##Define the optimization problem first
from pymoo.problems.functional import FunctionalProblem

problem=FunctionalProblem(n_var=2,objs=Perf,xl=np.array([1,10]),xu=np.array([5.5,200]))
problem.evaluate(np.array([3,150]))

##Hooke-Jeeves Search 
from pymoo.algorithms.soo.nonconvex.pattern import PatternSearch
from pymoo.optimize import minimize

resHK=minimize(problem,algorithm=PatternSearch(),verbose=True)
print("Best solution found: \nX = %s\nF = %s" % (resHK.X, resHK.F))


##Nelder-Mead Search 
from pymoo.algorithms.soo.nonconvex.nelder import NelderMead
resNM=minimize(problem,algorithm=NelderMead(),verbose=True)
print("Best solution found: \nX = %s\nF = %s" % (resNM.X, resNM.F))


(2)

#define a demand function
a=200
b=35

def Demand(Price,a=200,b=35):
  Dval=a-b*Price+stats.expon.rvs(scale=10,size=1)
  Dval=np.round(Dval)
  return max(Dval,0)

S=5000
salvage_value=0.5
shortage_cost=1
purchase_cost=1
def Perf(x):
  #resetting the seed is important!
  np.random.seed(5566)
  P=x[0]
  Q=x[1]
  profit=np.empty(S)
  for s in range(S):
    demand_s=Demand(P)
    profit[s]=P*min(Q,demand_s) - purchase_cost*Q + salvage_value*max(Q-demand_s,0) - shortage_cost*max(demand_s-Q,0)
  avgprofit=np.mean(profit)
  return -1*avgprofit
  
  ##Define the optimization problem first
from pymoo.problems.functional import FunctionalProblem

problem=FunctionalProblem(n_var=2,objs=Perf,xl=np.array([1,10]),xu=np.array([5.5,200]))
problem.evaluate(np.array([3,150]))

##Hooke-Jeeves Search 
from pymoo.algorithms.soo.nonconvex.pattern import PatternSearch
from pymoo.optimize import minimize

resHK=minimize(problem,algorithm=PatternSearch(),verbose=True)
print("Best solution found: \nX = %s\nF = %s" % (resHK.X, resHK.F))


##Nelder-Mead Search 
from pymoo.algorithms.soo.nonconvex.nelder import NelderMead
resNM=minimize(problem,algorithm=NelderMead(),verbose=True)


print("Best solution found: \nX = %s\nF = %s" % (resNM.X, resNM.F))


## Q2
Q2 (30%) 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, the demand depends on whether NCCU Griffins win or lose. If they win, we will continue to sell the shirt for \$25, and demand in the month afterwards is 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. All unsold shirts will be given away for free.

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

- Phase 1: sell "before" the Bowl game
- sale price: \$25
- demand distribution: Normal(μ=9000, σ=2000)
- Phase 2: sell "after" Bowl game
- win
- sale price: \$25
- demand distribution: Normal(μ=6000, σ=2000)
- lose
- sale price: \$12.5
- demand distribution: gamma distribution with mean 2000 and standard deviation of 1000

- salvage: \$0 (for free)

import numpy as np
import matplotlib.pyplot as plt
import random

# shirts cost&price
cost = 10
before_price = 25
win_price = 25
lose_price = 12.5
salvage = 0

# demend
random.seed(1206)
s = 10000
scale = (1000^2)/2000
shape = 2000/scale

shirts_demand_before = np.random.normal(9000, 2000, s)
shirts_demand_win = np.random.normal(6000, 2000, s)
shirts_demand_lose = np.random.gamma(shape, scale, s)

# win probability
win_prob = 0.4
win_or_not = np.random.choice([1,0],s,replace=True,p=[win_prob,1-win_prob])

# demand&profit 
expect_profit = np.empty(len(range(0,30500,500)))
expect_CVARq10 = np.empty(len(range(0,30500,500)))

for i in range(len(range(0,30500,500))):
  quantity = list(range(0,30500,500))
  profit = np.empty(s)
  for j in range(s):
      soldQ = min(quantity[i], shirts_demand_before[j])
      profit_before=(before_price-cost)*soldQ
      remainQ=quantity[i]-soldQ
      if win_or_not[j] == 1:
          profit_after = (win_price-cost)*min(remainQ,shirts_demand_win[j])-(cost-salvage)*max(remainQ-shirts_demand_win[j],0)
      else:
          profit_after = (lose_price-cost)*min(remainQ,shirts_demand_lose[j])-(cost-salvage)*max(remainQ-shirts_demand_lose[j],0)
      profit[j] = profit_before + profit_after
  expect_profit[i] = np.mean(profit)
  expect_CVARq10[i] = np.mean([x for x in profit if x <= np.quantile(profit, 0.1)])
  
# expected profits
plt.figure(figsize=(8,5))
plt.plot(quantity,expect_profit,color='red', linestyle='-',lw=5, alpha=0.5,label='expect profit')
plt.legend(loc = 'upper right')
plt.show()


print("The optimal production quantity to maximize expected profit = ",quantity[np.argmax(expect_profit)],"\nExpected profit = ",round(max(expect_profit),2),sep="")

The optimal production quantity to maximize expected profit = 12500
Expected profit = 143851.24

The optimal production quantity to maximize expected profit is 12,500 and the expected profit is 143,851.

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

# CVAR(q=10%)
plt.figure(figsize=(8,5))
plt.plot(quantity,expect_CVARq10,color='red', linestyle='-',lw=5, alpha=0.5,label='CVAR(q=10%)')
plt.legend(loc = 'upper right')
plt.show()


print("The optimal production quantity to maximize CVAR(q=10%) = ",quantity[np.argmax(expect_CVARq10)],"\nExpected profit = ",round(max(expect_CVARq10),2),sep="")

The optimal production quantity to maximize CVAR(q=10%) = 6900
Expected profit = 102170.0

The optimal production quantity to maximize CVAR(q=10%) is 6,900 and the expected profit is 102,170.

(3) Assuming we want to maximize expected profit in (1), what would be the expected value of perfect information about whether our local heroes will win the Bowl or not?

# perfect information
# before, win
cu_before_win = before_price - cost
co_before_win = cost-salvage
frac_before_win = cu_before_win / ( cu_before_win + co_before_win ) #0.6

# lose
cu_lose = lose_price - cost
co_lose = cost-salvage
frac_lose = cu_lose / ( cu_lose + co_lose ) #0.2

# fraction 
optimal_quantity_before = np.quantile(shirts_demand_before, frac_before_win)
optimal_quantity_win = np.quantile(shirts_demand_win, frac_before_win)
optimal_quantity_lose= np.quantile(shirts_demand_lose, frac_lose) 

# expect profit of perfect information
expect_profit_perfect =  np.empty(s)

for i in range(s):
  soldQ = min(shirts_demand_before[i],optimal_quantity_before)
  profit_before=(before_price-cost)*soldQ
  if win_or_not[i] == 1:
      profit_after = (win_price-cost)*min(shirts_demand_win[i],optimal_quantity_win)-(cost-salvage)*max(optimal_quantity_win-shirts_demand_win[i],0)
  else:
      profit_after = (lose_price-cost)*min(shirts_demand_lose[i],optimal_quantity_lose)-(cost-salvage)*max(optimal_quantity_lose-shirts_demand_lose[i],0)
  expect_profit_perfect[i] = profit_before + profit_after
  
print("Expected value of perfect information: ",round(np.mean(expect_profit_perfect),2),sep="")
print("maximum expected profit: ",round(max(expect_profit),2),sep="")
print("Benefit from perfect information: ",round(np.mean(expect_profit_perfect)-max(expect_profit),2),sep="")

* Expected value of perfect information: 156758.36
* maximum expected profit: 143851.24
* Benefit from perfect information: 12907.13

## Q3
Consider 𝑄7 with a possible range in [11, 110]. Please analyze strategies (a), (b), (c) below to find 𝑄7 Susan needs to have in order to maximize annual expected return E.

- Case Assumptions:
- Model of Demand for Analysts
Hi= Historic average analyst demand in month i
X= percentage unanticipated economic growth per year
X~ Normal(0%,5%)
Yi= Random noise in demand in month i compared to expected historic level in month i
Yi~ Normal(0%, 10%)
Di= Demand for analysts in month i
Di= Hi (1 + X) (1 + Yi)

- Model of Supply of Analysts
Qi= Number of offers made by bank to analysts to start in month i
Ai= Number of analysts who accept offer to start in month i
Ai~ Binomial( Qi, 0.7)
Ri= Percentage retention rate of analysts in month i
Ri~ Uniform(80%, 100%) for September and January
Ri~ Uniform(95%, 100%) for May, June, July and August
Ri~ Uniform(90%, 100%) for February, March, April, October, November and December
Pi= Number of analysts employed at start of month i
Pi+1= (Pi)(Ri) + Ai+1
P0= 63
- Monthly salary and indirect labor support cost
Monthly Salary $4,000
Monthly Indirect Cost $2,000
Monthly Total Cost $6,000

- Monthly Contribution to Earnings of productive analyst
Revenue per analyst month $10,000
Monthly Total Cost $6,000
Contribution per analyst $4,000

- Costs and Contribution of Transferred Workers
Inefficiencies of Transfers 60%
Monthly Total Cost of Transfer $9,600 (= $6,000 x 1.60)
Revenue per month $10,000
Monthly Total Cost $9,600
Contribution per transferfee $400

- Calculations of Earnings Contributions
Let Ei = contribution to earnings from analysts' work in month i
If demand equals supply of analysts then Ei = 4,000 Di = 4,000 Pi
But if not then earnings contributions are computed one of two ways each month, depending on whether there is a shortage or an excess of analysts in month i:
Excess of analysts in month i
*If Pi > Di then Ei = 10,000 Di - 6,000 Pi*
Shortage of analysts in month i
*If Di > Pi then Ei = (10,000 - 6,000) Pi + (10,000 - 9,600) (Di - Pi)*

(a) Take a Fixed-start strategy, i.e., send out Q7fixed offers to graduates in June, who will all be on board on July 1. Q7fixed? E7fixed?

import numpy as np
import matplotlib.pyplot as plt

#Demand simulation
S=10000
#historical demand started in April
H=[105, 95, 75, 70, 70, 110, 105, 90, 65, 80, 90, 120]
#percentage unanticipated economic growth per year
X= np.random.normal(0, 0.05, S)
#Random noise in demand in month i compared to expected historic level in month i
Y= np.random.normal(0, 0.1, S) 
#Di= Demand for analysts in month i
D= np.empty([12, S])
for i in range(12):
  D[i]= H[i]*(1 + X)*(1 + Y)
  
#Supply simulation
#fixed, started in April
#Ri= Percentage retention rate of analysts in month i 
R= np.empty([12, S])
for i in range(0, 12):
  if i==5|i==9:
    R[i]= np.random.uniform(0.8, 1, S) #for September and January
  elif(i==1|i==2|i==3|i==4): 
    R[i]= np.random.uniform(0.95, 1, S) #for May, June, July and August
  else:
    R[i]= np.random.uniform(0.9, 1, S) #for February, March, April, October, November and December
    
# objective function

#Pi= Number of analysts employed at start of month i
P= np.empty([12, S])
P[0,:]=63

def objective(Q_7):
  A_7= np.random.binomial(n= Q_7,p= 0.7,size= S)
  for j in range(1, 12):
    if j==3:
      P[j,]=P[(j-1),]*R[(j-1),]+A_7
    else:
      P[j,]=P[(j-1),]*R[(j-1),]
  E=np.empty([12, S])
  Exp= np.empty(12)
  ##
  for i in range(12):
    for j in range(S):
      if (P[i,j]>D[i,j]):
        E[i,j]= 10000*D[i,j]-6000*P[i,j]
      else:
        E[i,j]= 4000*P[i,j]+400*(D[i,j]-P[i,j])
    Exp[i]= np.mean(E[i,:])
  return np.round(np.sum(Exp),3)
  
Exp_return7=0
for Q in range(11, 111):
  if objective(Q)>Exp_return7:
    Exp_return7=objective(Q)
    Q_7=Q
print("Q7fixed=",Q_7,", E7fixed=", Exp_return7)

#Q7fixed= 57 , E7fixed= 3061422.828

(b) Take a Flexible-start strategy, i.e., send out 𝑄7flexible offers to graduates in June. Half will be on board on July 1, another half will be on board on September 1. 𝑄7flexible? 𝐸7flexible?

# objective function #flexible
S=10000

#Pi= Number of analysts employed at start of month i
P= np.empty([12, S])
P[0,:]=63

def objective_f(Q_7):
  A_7= np.random.binomial(n= Q_7,p= 0.7,size= S)
  for j in range(1, 12):
    if (j==3 or j==5): #July#September
      P[j,]=P[(j-1),]*R[(j-1),]+0.5*A_7
    else:
      P[j,]=P[(j-1),]*R[(j-1),]
  E=np.empty([12, S])
  Exp= np.empty(12)
  ##
  for i in range(12):
    for j in range(S):
      if (P[i,j]>D[i,j]):
        E[i,j]= 10000*D[i,j]-6000*P[i,j]
      else:
        E[i,j]= 4000*P[i,j]+400*(D[i,j]-P[i,j])
    Exp[i]= np.mean(E[i,:])
  return np.round(np.sum(Exp),3)
  
Exp_return7f=0
for Q in range(11, 111):
  if objective_f(Q)>=Exp_return7f:
    Exp_return7f=objective_f(Q)
    Q_7=Q
print("Q7flexible=",Q_7, ", E7flexible=", Exp_return7f)

#Q7flexible= 70 , E7flexible= 3307391.885

( c ) Add a recruiting event in December and send out 𝑄12 offers in [11,110], who will be on board on January 1. This event will operate with the Fixed-start strategy or the Flexible-start strategy. What will be the new 𝑄12fixed? 𝐸12fixed? How about 𝑄12flexible? 𝐸12flexible? Based on your simulation, will the December recruiting event increase expected return?

Ans:
Case中 "Susan also thought there might be some merit to engaging in December recruiting in order to capture students who graduate at mid-year. Susan thought that mid-year recruiting, if done correctly, could help to alleviate analyst staffing shortages. In fact, the option to do December recruiting might even allow Casterbridge to hire fewer analysts for July 1 and then make a mid-year correction if need be . Analysts hired in the December round would be asked to start work at the beginning of January."

但我們的模擬從4月開始(因為P0)，所以無法知道12月招募所帶來的效益
另外，因Case中未提及flexible的12月招募計畫應該從哪一個月開始，因此我們沿用上題討論在4月之前的另一半員工應該在哪個月進入

(1)fixed

#Demand simulation
S=10000
#historical demand started in April
H=[105, 95, 75, 70, 70, 110, 105, 90, 65, 80, 90, 120]
#percentage unanticipated economic growth per year
X= np.random.normal(0, 0.05, S)
#Random noise in demand in month i compared to expected historic level in month i
Y= np.random.normal(0, 0.1, S) 
#Di= Demand for analysts in month i
D= np.empty([12, S])
for i in range(12):
  D[i]= H[i]*(1 + X)*(1 + Y)

#Supply simulation
#fixed, started in April
#Ri= Percentage retention rate of analysts in month i 
R= np.empty([12, S])
for i in range(0, 12):
  if i==5|i==9:
    R[i]= np.random.uniform(0.8, 1, S) #for September and January
  elif(i==1|i==2|i==3|i==4): 
    R[i]= np.random.uniform(0.95, 1, S) #for May, June, July and August
  else:
    R[i]= np.random.uniform(0.9, 1, S) #for February, March, April, October, November and December

# objective function

#Pi= Number of analysts employed at start of month i
P= np.empty([12, S])
P[0,:]=63

def objective(Q_12):
  A_7= np.random.binomial(n= 57,p= 0.7,size= S) #使用(a)答案模擬 #57
  A_12= np.random.binomial(n= Q_12,p= 0.7,size= S)
  for j in range(1, 12):
    if j==3:
      P[j,]=P[(j-1),]*R[(j-1),]+A_7
    elif j==9: #Jan
      P[j,]=P[(j-1),]*R[(j-1),]+A_12
    else:
      P[j,]=P[(j-1),]*R[(j-1),]
  E=np.empty([12, S])
  Exp= np.empty(12)
  ##
  for i in range(12):
    for j in range(S):
      if (P[i,j]>D[i,j]):
        E[i,j]= 10000*D[i,j]-6000*P[i,j]
      else:
        E[i,j]= 4000*P[i,j]+400*(D[i,j]-P[i,j])
    Exp[i]= np.mean(E[i,:])
  return np.round(np.sum(Exp),3)

Exp_return12=0
for Q in range(11, 111):
  if objective(Q)>Exp_return12:
    Exp_return12=objective(Q)
    Q_12=Q
print("Q12fixed=",Q_12,", E12fixed=", Exp_return12)

#Q12fixed= 28 , E12fixed= 3142560.02

(2) flexible
假設招募計畫的另一半人入職時間從2~4月皆有可能，研究哪個月入職最好

# objective function

#Pi= Number of analysts employed at start of month i
P= np.empty([12, S])
P[0,:]=63

def objective_f(Q_12):
  A_7= np.random.binomial(n= 57,p= 0.7,size= S) #使用(a)答案模擬 #57
  A_12= np.random.binomial(n= Q_12,p= 0.7,size= S)
  for j in range(1, 12):
    if j==3:
      P[j,]=P[(j-1),]*R[(j-1),]+A_7
    elif j==9 or j== 11: #Jan #Mar
      P[j,]=P[(j-1),]*R[(j-1),]+0.5*A_12
    else:
      P[j,]=P[(j-1),]*R[(j-1),]
  E=np.empty([12, S])
  Exp= np.empty(12)
  ##
  for i in range(12):
    for j in range(S):
      if (P[i,j]>D[i,j]):
        E[i,j]= 10000*D[i,j]-6000*P[i,j]
      else:
        E[i,j]= 4000*P[i,j]+400*(D[i,j]-P[i,j])
    Exp[i]= np.mean(E[i,:])
  return np.round(np.sum(Exp),3)

Exp_return12_f=0
for Q in range(11, 111):
  if objective(Q)>Exp_return12_f:
    Exp_return12_f=objective(Q)
    Q_12=Q

print("Q12flexible=",Q_12,"Month=3, E12flexible=", Exp_return12_f)

#Q12flexible= 28 Month=2, E12flexible= 3142099.742
#Q12flexible= 28 Month=3, E12flexible= 3143662.554
另一半的員工在3月時加入

res = {"E7fixed":Exp_return7,
    'E7flexible':Exp_return7f,
    'E12fixed':Exp_return12,
    'E12flexible':Exp_return12_f,
    'Max':max(Exp_return7,Exp_return7f,Exp_return12,Exp_return12_f)}
print(res)
print("Expected return最大的方案是7月招募的flexible計畫，Expected return=",max(Exp_return7,Exp_return7f,Exp_return12,Exp_return12_f))

{'E7fixed': 3061422.828,
'E7flexible': 3307391.885,
'E12fixed': 3142560.02,
'E12flexible': 3143662.554,
'Max': 3307391.885}

Expected return最大的方案是7月招募的flexible計畫，Expected return= 3307391.885