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

#  Required packages
require(statmod)
require(randtoolbox)

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

####################################################
#   Chunk 1: Use Gauss-Hermite quadrature and
#            check moments, up to kurtosis
mu <- 50
sigma <- 20
N <- 5
out <- gauss.quad(N, kind="hermite")
w1 <- out$weights/sqrt(pi)
x1 <- sqrt(2)*sigma*out$nodes + mu
out <- gauss.quad.prob(N,dist="normal",mu=mu,sigma=sigma)
w2 <- out$weights
x2 <- out$nodes 
rbind(w1,w2)
rbind(x1,x2)

sum(w1)
sum(w1*x1)
sqrt(sum(w1*x1^2)-sum(w1*x1)^2)
sum(w1*(x1-50)^3)/20^3
sum(w1*(x1-50)^4)/20^4


####################################################
#   Chunk 2: Use Gauss-Hermite quadrature to find 
#            the expected value of exp(X), with X normal
mu <- 4
sigma <- 2
numNodes <- c(5,10,15,20)
trueValue <- exp(mu + sigma^2/2)
for (i in 1:length(numNodes)){
  out <- gauss.quad(numNodes[i], kind="hermite")
  w <- out$weights/sqrt(pi)
  x <- sqrt(2)*sigma*out$nodes + mu
  approxValue <- sum(w*exp(x))
  percError <- 100*abs(trueValue-approxValue)/trueValue
  cat("N=",numNodes[i]," True=",trueValue, " Approx=",
       approxValue, " %error", percError,"\n")
}
 
####################################################
#   Chunk 3: Naive Monte Carlo integration
#
###   Show variability  
set.seed(55555)
mean(exp(runif(10)))
mean(exp(runif(10)))
mean(exp(runif(10)))
mean(exp(runif(100)))
mean(exp(runif(100)))
mean(exp(runif(100)))
mean(exp(runif(1000000)))

###   Compare against Gaussian quadrature 
out=gauss.quad.prob(10,dist="uniform",l=0,u=1)
sum(out$weight*exp(out$nodes))

###   Find confidence intervals
set.seed(55555)
out <- t.test(exp(runif(100)))
out$estimate
out$conf.int

set.seed(55555)
X=exp(rnorm(100,mean=4,sd=2))
t.test(X)
  
####################################################
#   Chunk 4: Compare random and Halton point sets
#
set.seed(55555)
plot(runif(100),runif(100),pch=20)
grid(nx=10)
x <- halton(100,dim=2)
plot(x[,1], x[,2], pch=20)
grid(nx=10)

####################################################
#   Chunk 5: Display Korobov point sets

### Case 1
a <- 3
d <- 2
N <- 10
X <- korobov(a,d,N)
X
plot(X[,1],X[,2],pch=20)
sort(X[,1])
sort(X[,2])

### Case 2
N <- 12
X <- korobov(a,d,N)
sort(X[,1])
sort(X[,2])

####################################################
#   Chunk 6: R library functions for numerical quadrature

### Straighforward case
f <- function(x) exp(-x)*sin(10*x)
integrate(f, lower=0, upper=2*pi)

### Illustrate vectorization issues
g <- function(x) return(1)
integrate(g,0,10)   # WRONG example  !!!!
g1 <- function(x) rep(1,length(x))
g(1:5)
g1(1:5)
integrate(g1,0,10)
g2 <- Vectorize(g)
g2(1:5)
integrate(g2,0,10)

### Issues with infinite intervals
integrate(dnorm,0,2)
integrate(dnorm,0,2000)
integrate(dnorm,0,20000000)
integrate(dnorm,0,Inf)

### A multidimensional example
library(cubature)
h <- function(x) exp(-x[1]*x[2])*(sin(6*pi*x[1])+
                                   cos(8*pi*x[2]))
adaptIntegrate(h,lower=c(0,0),upper=c(1,1))$integral

