####################################################
#    Scripts for Chapter 10 of
#    P. Brandimarte
#    Handbook in Monte Carlo Simulation
#    Applications in Financial Engineering,
#    Risk Management, and Economics
#    Wiley, 2014

#  Required packages
require(polynom)
require(linprog)
require(quadprog)
require(GenSA)
require(GA)
require(pso)
require(statmod)
require(randtoolbox)

#  Source the function file
source("Chapter10_Functs.R")

####################################################
#   Chunk 1: The role of convexity
####################################################

### Optimizing a polynomial may result in local minima
p <- c(30, -59.5, 39, -10.5, 1)
f <- as.function(polynomial(p))
x <- seq(from=1,to=4,by=0.01)
plot(x,f(x),type='l')
out1 <- optim(par=4,fn=f,method="BFGS")
out1$par; out1$value
out2 <- optim(par=1,fn=f,method="BFGS")
out2$par; out2$value
xx <- seq(from=1,to=4,by=0.01)
plot(xx,f(xx),type='l')

####################################################
#   Chunk 2: Linear and quadratic programming
####################################################

### Example 1: Using linprog
Amat <- rbind(c(15,35),c(25,15),c(1,0),c(0,1))
bvet <- c(2400,2400,100,50)
cvet <- c(45,60)
result <- solveLP(cvet,bvet,Amat,maximum=TRUE)
cat('optimal obj ', result$opt, '\n')
cat('optimal x ', result$solution, '\n')

### Example 2: Using quadprog
Qmat <- rbind(c(2,0),c(0,2))
bvet <- c(4)
Amat <- matrix(c(1,1),nrow=2)
result <- solve.QP(Qmat, c(0,0),Amat,bvet,meq=1)
cat('optimal obj ', result$value, '\n')
cat('optimal x ', result$solution, '\n')


####################################################
#   Chunk 3: Simulated annealing
####################################################

### Example 1: Compare classical nonlinear programming
###            and homemade simulated annealing
Rastrigin<-function(x){sum(x^2-10*cos(2*pi*x))+10*length(x)}
dim <- 20
lb <- rep(-100,dim)
ub <- rep(100,dim)
set.seed(55555)
x0 <- runif(dim,min=lb,max=ub)
out <- optim(par=x0, fn=Rastrigin, lower=lb, upper=ub, 
                       method="L-BFGS-B")
out[c("value","counts")]
set.seed(55555)
out <- NaiveSA(fn=Rastrigin,par=x0,temp0=100,
                    beta=0.9999,maxit=200000,sigma=1)
out[c("value","counts")]

### Example 2: Compare better implementation of SA
###            and generalized simulated annealing
set.seed(55555)
out <- optim(par=x0,fn=Rastrigin,method="SANN",
             control=list(temp=100,tmax=10,maxit=200000))
out[c("value","counts")]
set.seed(55555)
out <- GenSA(fn = Rastrigin,lower=lb,upper=ub,par=x0)
out[c("value","par","counts")]


####################################################
#   Chunk 4: Using genetic algorithms
####################################################

Rastrigin<-function(x){sum(x^2-10*cos(2*pi*x))+10*length(x)}
dim <- 20
lb <- rep(-100,dim)
ub <- rep(100,dim)
Fitness <- function(x) -Rastrigin(x)
set.seed(55555)
out <- ga(type="real-valued",fitness=Fitness,min=lb,max=ub,
          popSize=2000,maxiter=200)
summary(out)


####################################################
#   Chunk 5: Using particle swarm optimization
####################################################

Rastrigin<-function(x){sum(x^2-10*cos(2*pi*x))+10*length(x)}
dim <- 20
lb <- rep(-100,dim)
ub <- rep(100,dim)
set.seed(55555)
x0 <- runif(dim,min=lb,max=ub)
out <- psoptim(par=x0,fn=Rastrigin,lower=lb,upper=ub,
               control=list(maxit=10000,trace=1,REPORT=1000,
                            hybrid=F))
out <- psoptim(par=x0,fn=Rastrigin,lower=lb,upper=ub,
               control=list(maxit=1000,trace=1,REPORT=100,
                            hybrid=T))

####################################################
#   Chunk 6: Illustrating scenario generation and 
#            stability issues in stochastic optimization
#            by a simple newsvendor model
####################################################

### Set parameters and solve optimally
pur <- 10
sell <- 12
mdown <- 5
mu <- 100
sigma <- 30
newsOptimal(pur,sell,mdown,mu,sigma)

### Numerical solution by stochastic LP
### use naive MC and check in-sample stability
m <- sell-pur
cu <- pur-mdown
set.seed(55555)
numScen <- 50
probs <- rep(1/numScen,numScen)
qv <- numeric(5)
for (k in 1:5){
  demand <- pmax(rnorm(numScen,mu,sigma),0)
  out <- newsLP(m,cu,demand,probs)
  qv[k] <- out$q
  cat('obj ',round(out$cost,4),' q ',round(out$q,1),'\n')
}
# Check out-of-sample stability
outScost <- function(q,m,cu,demand){
  mean(cu*pmax(q-demand,0)+m*pmax(demand-q,0))}
set.seed(55555)
outsDemand <- pmax(rnorm(100000,mu,sigma),0)
for (k in 1:5){
  trueCost <- outScost(qv[k],m,cu,outsDemand)
  cat('true cost ',round(trueCost,4),'\n')
}

### Use Gaussian quadrature
quad <- gauss.quad.prob(10,dist="normal",mu=mu,sigma=sigma)
out <- newsLP(m,cu,pmax(quad$nodes,0),quad$weights)
cat('obj ',round(out$cost,4),' q ',round(out$q,1),'\n')
trueCost <- outScost(out$q,m,cu,outsDemand)
cat('true cost ',round(trueCost,4),'\n')
quad$nodes
quad$weights

### Use stratification (sort of...)
numVals <- 50
aux <- (0:numVals)/numVals
U <- (aux[-1]+aux[-(numVals+1)])/2
demand <- pmax(0, mu+sigma*qnorm(U))
probs <- rep(1/numVals,numVals)
out <- newsLP(m,cu,demand,probs)
cat('obj ',round(out$cost,4),' q ',round(out$q,1),'\n')
trueCost <- outScost(out$q,m,cu,outsDemand)
cat('true cost ',round(trueCost,4),'\n')
range(demand)

### Use Sobol sequences
numVals <- 50
eps <- sobol(numVals,dim=1,normal=TRUE,scrambling=2)
demand <- pmax(0, mu+sigma*eps)
probs <- rep(1/numVals,numVals)
out <- newsLP(m,cu,demand,probs)
cat('obj ',round(out$cost,4),' q ',round(out$q,1),'\n')
trueCost <- outScost(out$q,m,cu,outsDemand)
cat('true cost ',round(trueCost,4),'\n')


####################################################
#   Chunk 7: Illustrating numerical dynamic programming
#            and the approximation of value functions 
#            by cubic splines in a determinstic 
#            resource allocation problem
####################################################

### Illustration of cubic splines
x <- 1:5
y = c(8,5,3,7,4)
f <- splinefun(x,y)
plot(x,y)
pts <- seq(from=1,to=5,by=0.01)
lines(pts,f(pts))

### Solving the problem 
budget <- 20
numPoints <- 50
f <- function(x) sqrt(x)
funcList <- list(f,f,f)
out <- Budgeting(budget,funcList,numPoints)
out$par
out$value


####################################################
#   Chunk 8: Illustrating value iteration for a
#            very simple MDP
####################################################

### Case 1
probs = c(.7,.8,.9,0)
payoff = c(0,10,20,30)
discount=.8
out <- PW311_ValueIteration(probs, payoff, discount)
out

### Case 2
probs = c(.6,.6,.6,0)
payoff = c(0,10,20,30)
discount=.95
out <- PW311_ValueIteration(probs, payoff, discount)
out

### Case 3
probs = c(.6,.6,.6,0)
payoff = c(0,10,20,30)
discount=.99
out <- PW311_ValueIteration(probs, payoff, discount)
out

####################################################
#   Chunk 9: Solving the consumption-saving problem
#            by numerical stochastic DP
####################################################

### Experiment 1: Value functions for logarithmic utility
utilFun <- function(x) log(x); gamma <- 0.97
numSteps <- 10; W0 <- 100
mu <- 0.07; sigma <- 0.2; rf <-0.03
empl0 <- 2; income <- c(5,20,40)
transMatrix <- matrix(
  c(0.6,0.3,0.1,0.2,0.5,0.3,0.2,0.1,0.7),
  nrow=3,byrow=TRUE)
Wmin <- 1; Wmax <- 500; numPoints <- 30
numScenarios <- 11
scen <- MakeScenariosGQ(mu,sigma,numScenarios)
retScenarios <- scen$retScenarios
probs <- scen$probs
valFuns <- DPSaving(W0,empl0,utilFun,gamma,numSteps,
                    numPoints,Wmin,Wmax,income,transMatrix,
                    retScenarios,probs,rf)

x=Wmin:Wmax
par(mfrow=c(2,2))
plot(x,valFuns[[10]][[2]](x),type='l',ylab='time=10')
plot(x,valFuns[[9]][[2]](x),type='l',ylab='time=9')
plot(x,valFuns[[5]][[2]](x),type='l',ylab='time=5')
plot(x,valFuns[[1]][[2]](x),type='l',ylab='time=1')
par(mfrow=c(1,1))


### Experiment 2: Value functions for power utility
###       and different splines (code NOT shown in book)
eta <- 5
utilFun2 <- function(x) x^(1-eta)/(1-eta)
# Natural splines
valFunsNAT <- DPSaving(W0,empl0,utilFun2,gamma,numSteps,numPoints,
                       Wmin,Wmax,income,transMatrix,retScenarios,
                       probs,rf)
# Monotonicity preserving splines
valFuns2 <- DPSaving(W0,empl0,utilFun2,gamma,numSteps,numPoints,
                     Wmin,Wmax,income,transMatrix,retScenarios,
                     probs,rf,"monoH.F")
x=Wmin:(Wmax/3)
par(mfrow=c(2,2))
plot(x,valFunsNAT[[10]][[2]](x),type='l',ylab='time=10,NAT')
plot(x,valFunsNAT[[9]][[2]](x),type='l',ylab='time=9,NAT')
plot(x,valFuns2[[10]][[2]](x),type='l',ylab='time=10,MH')
plot(x,valFuns2[[9]][[2]](x),type='l',ylab='time=9,MH')
par(mfrow=c(1,1))


### Experiment 3: Comparing a fixed and a
### DP policy (out of sample, log utility)
numRepl <- 100
set.seed(55555)
scenarios <- MakeSamplePaths(empl0,transMatrix,mu,sigma,
                             numSteps,numRepl)
alpha <- 0.5
fractCons <- 0.5
out <- SimulateFixedSaving(W0,empl0,utilFun,gamma,
                           income,transMatrix,mu,sigma,rf,
                           scenarios,alpha, fractCons)
utilValsFixed <- out$utilVals
out <- SimulateDPSaving(W0,empl0,utilFun,gamma,valFuns,
                        income,transMatrix,retScenarios,
                        probs,rf,scenarios)
utilValsDP <- out$utilVals
mean(utilValsFixed)
mean(utilValsDP)
t.test(utilValsDP,utilValsFixed)

### Experiment 4: Comparing a fixed and a
### DP policy (out of sample, power utility)
### Code NOT shown in book. This experiment gives slightly 
### different results with respect to those shown in the book, 
### which were actually obtained using 10 scenarios in GQ
numRepl <- 100
set.seed(55555)
scenarios <- MakeSamplePaths(empl0,transMatrix,mu,sigma,
                             numSteps,numRepl)
alpha <- 0.5
fractCons <- 0.5
out <- SimulateFixedSaving(W0,empl0,utilFun2,gamma,
                           income,transMatrix,mu,sigma,rf,
                           scenarios,alpha,fractCons)
utilValsFixed <- out$utilVals
out <- SimulateDPSaving(W0,empl0,utilFun2,gamma,valFuns2,
                        income,transMatrix,retScenarios,
                        probs,rf,scenarios)
utilValsDP <- out$utilVals
mean(utilValsFixed)
mean(utilValsDP)
t.test(utilValsDP,utilValsFixed)


####################################################
#   Chunk 10: Solving a simple MDP by Q-learning
####################################################

### Experiment 1: compare Q-learning and value iteration
probs <- c(.7,.8,.9,0); payoff <- c(0,10,20,30)
discount <- 0.8
PW311_ValueIteration(probs, payoff, discount)
numSteps <- 10000; startState <- 1; alpha <- 0.1
startQFactors <- matrix(c(rep(50,4),payoff),ncol=2)
set.seed(55555)
PW311_QLearning(probs,payoff,discount,numSteps,
                  startState,startQFactors,alpha)

### Experiment 2: Change the initialization of Q-factors
###               and see slow convergence
startQFactors <- matrix(c(rep(500,4),payoff),ncol=2)
set.seed(55555)
PW311_QLearning(probs,payoff,discount,numSteps,
                   startState,startQFactors,alpha)

### Experiment 3: Change the initialization of Q-factors
###               and see lack of convergence
startQFactors <- matrix(c(rep(1,4),rep(0,4)),ncol=2)
set.seed(55555)
PW311_QLearning(probs,payoff,discount,numSteps,
                   startState,startQFactors,alpha)

### Experiment 4: Improve exploration by the eps-greedy strategy
startQFactors <- matrix(c(rep(1,4),rep(0,4)),ncol=2)
eps <- 0.1
set.seed(55555)
PW311_QepsLearning(probs,payoff,discount,numSteps,
                     startState,startQFactors,alpha,eps)
