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

#  Required packages
require(MASS) # this package includes mvrnorm 

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

####################################################
#   Chunk 1: Shortfall probability for a portfolio
#            of vanilla calls and puts
####################################################

### Assessing the accuracy of delta-gamma approximation
S0 <- 100; r <- 0.05; mu <- 0.05
sigma <- 0.3; T <- 0.5; K <- 100
numStocks <- 10
numCalls <- -10; numPuts <- -5
numRepl <- 3000
dt <- 0.1
thresh <- 50
set.seed(55555)
out1 <- ProbShortfallCrudeMC(S0,K,T,mu,r,sigma,numStocks,
                             numCalls,numPuts,dt,numRepl,thresh)
set.seed(55555)
out2 <- ProbShortfallDG(S0,K,T,mu,r,sigma,numStocks,
                        numCalls,numPuts,dt,numRepl,thresh)
out1$estimate
out2$estimate
plot(out1$loss,out2$loss,pch=20)
cor(out1$loss,out2$loss)
cor(out1$loss>50,out2$loss>50)
cor(out1$loss>150,out2$loss>150)
cor(out1$loss>250,out2$loss>250)


### Assessing the reduction in variance by exponential tilting
S0 <- 100; r <- 0.05; mu <- 0.05; sigma <- 0.3
T <- 0.5; K <- 100; numStocks <- 10
numCalls <- -10; numPuts <- -5
numRepl <- 3000
dt <- 0.1; thresh <- 250
set.seed(55555)
N <- 100
outMC <- numeric(N)
outIS <- numeric(N)
for (k in 1:N){
  aux <- ProbShortfallCrudeMC(S0,K,T,mu,r,sigma,
                              numStocks,numCalls,numPuts,dt,numRepl,thresh)
  outMC[k] <- aux$estimate
  outIS[k] <- ProbShortfallIS(S0,K,T,mu,r,sigma,numStocks,
                              numCalls,numPuts,dt,numRepl,thresh)
}
# get estimates and confidence intervals
estMC <- t.test(outMC)
estIS <- t.test(outIS)
percMC<-(estMC$conf.int[2]-estMC$conf.int[1])/estMC$estimate
percIS<-(estIS$conf.int[2]-estIS$conf.int[1])/estIS$estimate
cat("estimate MC =",estMC$estimate,"\n")
cat("estimate IS =",estIS$estimate,"\n")
cat("percentage error crude MC   =",percMC*100,"%\n")
cat("percentage error imp. samp. =",percIS*100,"%\n")


####################################################
#   Chunk 2: Plotting V@R as a function of portfolio 
#   weights: The case of two assets and discretized scenarios
####################################################

s1 <- 0.4
s2 <- 0.3
rho <- 0.6
Sigma <- matrix(c(s1^2,rho*s1*s2,rho*s1*s2,s2^2), nrow=2)
mu <- c(0.12,0.10)

### define the auxiliary function here
plotVaR <- function(numScen,alpha){
  rets <- mvrnorm(n=numScen,mu=mu,Sigma=Sigma)
  w <- seq(from=0,to=1,by=0.001)
  n <- length(w)
  VaR <- numeric(n)
  for (i in 1:n){
    portRets <- w[i]*rets[,1]+(1-w[i])*rets[,2]
    sortedRets <- sort(portRets)
    # use naive estimation of quantiles, based on
    # sorted ordered statistics
    VaR[i] <- -sortedRets[alpha*numScen]
  }
  plot(w,VaR,type='l')
}

set.seed(55555)
par(mfrow=c(2,2))
plotVaR(100,0.05)
plotVaR(100,0.05)
plotVaR(100,0.05)
plotVaR(100,0.05)
par(mfrow=c(1,1))
plotVaR(10000,0.05)


####################################################
#   Chunk 3: Estimating the cost of hedging strategies
#   for a vanilla call option 
####################################################

### Experiment 1: Stop-loss strategy, and checking
### the effect of code vectorization
S0 <- 50; K <- 52; mu <- 0.1;
sigma <- 0.4; r <- 0.05; T <- 5/12;
numRepl <- 100000; numSteps <- 10
set.seed(55555)
paths <- simvGBM(S0,mu,sigma,T,numSteps,numRepl)
system.time(cost <- stopLoss(paths,K,r,T))
system.time(costV <- stopLossV(paths,K,r,T))
cost$mean
costV$mean

### Experiment 2: Delta-hedging strategy vs. stop-loss
S0 <- 50; K <- 52; mu <- 0.1;
sigma <- 0.4; r <- 0.05; T <- 5/12;
numRepl <- 100000; numSteps <- 10
# exact call price
C <- EuVanillaCall(S0,K,T,r,sigma)
# try 10 steps
set.seed(55555)
paths <- simvGBM(S0,mu,sigma,T,numSteps,numRepl)
SL <- stopLossV(paths,K,r,T)
DC <- deltaHedging(paths,K,sigma,r,T)
cat("true price = ", C, "\n")
cat("StopLoss 10 steps: mean =", SL$mean,
    "sd =", SL$sd, "\n")
cat("DeltaHedging 10 steps: mean =", DC$mean,
    "sd =", DC$sd, "\n")
# try 100 steps
numSteps <- 100
set.seed(55555)
paths <- simvGBM(S0,mu,sigma,T,numSteps,numRepl)
SL <- stopLossV(paths,K,r,T)
DC <- deltaHedging(paths,K,sigma,r,T)
cat("StopLoss 100 steps: mean =", SL$mean,
    "sd =", SL$sd, "\n")
cat("DeltaHedging 100 steps: mean =", DC$mean,
    "sd =", DC$sd, "\n")

