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

#  Required packages
require(randtoolbox)
require(cubature)
require(bitops)

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

####################################################
#   Chunk 1: Halton sequences
####################################################

### Example 1: Illustration of Halton sequences
halton(10, dim = 1)
halton(10, dim = 3)

### Example 2: Comparing Halton and pseudorandom seqs
set.seed(55555)
haltonSeq <- halton(1000, dim = 2)
windows()
plot(haltonSeq[,1],haltonSeq[,2],pch=20)
grid(nx=10,col='black')
windows()
plot(runif(1000),runif(1000),pch=20)
grid(nx=10,col='black')

### Example 3: Numerical integration
# using classical quadrature
f <- function(x) exp(-x[1]*x[2])*
  (sin(6*pi*x[1])+cos(8*pi*x[2]))
adaptIntegrate(f,rep(0,2),rep(1,2),tol=1e-4)
# using pseudorandom sampling
set.seed(55555)
mean(apply(matrix(runif(20000),ncol=2),1,f))
mean(apply(matrix(runif(20000),ncol=2),1,f))
mean(apply(matrix(runif(20000),ncol=2),1,f))
mean(apply(matrix(runif(20000),ncol=2),1,f))
mean(apply(matrix(runif(20000),ncol=2),1,f))
# using Halton seqs
mean(apply(halton(10000,dim=2),1,f))

### Example 4: Issues with high-dimensional Halton seqs
X <- halton(n=1000,dim=30)
plot(X[,30],pch=20)
1/X[1,30]
plot(X[,30],X[,29],pch=20)


####################################################
#   Chunk 2: Sobol sequences
####################################################

### Example 1: Illustrating Sobol sequences
sobol(n=10,dim=4)

### Example 2: Illustrating bit-level operations
###            and generation of direction numbers
m <- c(1, 3, 7)
i <- 4
m[i] <- bitXor(4*m[i-2],bitXor(8*m[i-3],m[i-3]))
i <- 5
m[i] <- bitXor(4*m[i-2],bitXor(8*m[i-3],m[i-3]))
i <- 6
m[i] <- bitXor(4*m[i-2],bitXor(8*m[i-3],m[i-3]))
m
p <- c(1,0,1,1) # a primitive polynomial
m0 <- c(1,3,7)
GetDirNumbers(p,m0,6)

### Example 3: Illustrating Gray codes
gray <- function (x) bitXor(x,bitShiftR(x,b=1))
codes <- matrix(0,nrow=16,ncol=4)
for (i in 1:16) codes[i,] <- 
     as.numeric(intToBits(gray(i-1)))[c(4,3,2,1)]
codes


####################################################
#   Chunk 3: Scrambling low-discrepancy seqs
####################################################

### Example 1: A scrambled Sobol sequence
sobol(n=10,dim=4,scrambling=1)

### Example 2: The impact of scrambling on integration
# same f as above, repeated for conveniency
f <- function(x) exp(-x[1]*x[2])*
  (sin(6*pi*x[1])+cos(8*pi*x[2]))
adaptIntegrate(f,rep(0,2),rep(1,2),tol=1e-6)$integral
mean(apply(sobol(10000,dim=2),1,f))
mean(apply(sobol(10000,dim=2,scrambling=1),1,f))
mean(apply(sobol(10000,dim=2,scrambling=2),1,f))


####################################################
#   Chunk 4: Option pricing and path generation
####################################################

### Example 1: Pricing a vanilla call
S0 <- 50; K <- 50; sigma <- 0.4
r <- 0.05; T <- 5/12
EuVanillaCall(S0,K,T,r,sigma)
set.seed(55555)
BlsLowD(S0,K,T,r,sigma,runif(5000))
BlsLowD(S0,K,T,r,sigma,halton(5000,dim=1))
BlsLowD(S0,K,T,r,sigma,sobol(5000,dim=1))
BlsLowD(S0,K,T,r,sigma,
        sobol(5000,dim=1,scrambling=1))

### Example 2: Sampling a GBM process
S0 <- 50; mu <- 0.1; sigma <- 0.4;
T <- 1/12; numSteps <- 30; numRepl <- 5
pathsH <- lowDiscrGBM(S0,mu,sigma,T,numSteps,
                        numRepl,type="halton")
pathsS <- lowDiscrGBM(S0,mu,sigma,T,numSteps,numRepl)
par(mfrow=c(1,2))
plot(pathsH[1,],type='l',ylim=range(pathsH),
       xlab="Halton",ylab="")
for (k in 2:numRepl) lines(pathsH[k,])
plot(pathsS[1,],type='l',ylim=range(pathsS),
       xlab="Sobol",ylab="")
for (k in 2:numRepl) lines(pathsS[k,])
par(mfrow=c(1,1))

### Some more paths (code not included in the text)
numRepl <- 10
paths <- lowDiscrGBM(S0,mu,sigma,T,numSteps,numRepl,scrambling=1)
plot(paths[1,],type='l',ylim=range(paths),xlab="Scrambled Sobol",ylab="")
for (k in 2:numRepl) lines(paths[k,])

