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

#  Required packages
require(MASS)
require(mvtnorm)
require(copula)

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

####################################################
#   Chunk 1: Linear congruential generators
####################################################

### Running an illustrative example 
a <- 5
c <- 3
m <- 16
seed <- 7
N <- 20
seq <- LCG(a,c,m,seed,N)
# unpack the output list into two sequences
USeq <- seq$U 
ZSeq <- seq$Z
# sprintf generates an output string by a C-like format
# specification, which is concatenated with a newline 
# character and printed by cat
for (k in 1:N)
  cat(sprintf('%3d %3d %7.4f', k, ZSeq[k], USeq[k]), "\n")

### Illustration of the lattice structure of LCGs
m <- 2048
a <- 65
c <- 1
seed <- 0
U <- LCG(a,c,m,seed, 2048)$U
# first pair of plots
par(mfrow = c(2,1))
plot(U[1:(m-1)], U[2:m], pch=20)
plot(U[1:511], U[2:512], pch=20)
# second plot
a <- 1365
c <- 1
U <- LCG(a,c,m,seed, 2048)$U
win.graph()
plot(U[1:(m-1)], U[2:m], pch=20)
par(mfrow = c(1,1))


####################################################
#   Chunk 2: Alternative generators
####################################################

### How to check the internal state of a 
### Mersenne twister generator
mode(.Random.seed)
class(.Random.seed)
length(.Random.seed)
set.seed(55555)
cat(.Random.seed[1:6],"\n")
runif(3)
cat(.Random.seed[1:6],"\n")
set.seed(55555)
cat(.Random.seed[1:6],"\n")

### Changing to a Wichmann--Hill generator
RNGkind(kind="Mersenne-Twister")
set.seed(55555)
cat(.Random.seed[1:6],"\n")
RNGkind(kind="Wichmann-Hill")
set.seed(55555)
cat(.Random.seed,"\n")
RNGkind(kind="default")
set.seed(55555)
cat(.Random.seed[1:6],"\n")

####################################################
#   Chunk 3: The inverse transform method
####################################################

### Sampling from an empirical discrete distribution
set.seed(55555)
values <- 1:5
probs <- c(0.1, 0.3, 0.4, 0.15, 0.05)
sample <- EmpiricalDrnd(values,probs,10000)
freqs <- hist(sample,breaks=0:5,plot=FALSE)
plot(freqs$counts,type='h',lwd=15)


####################################################
#   Chunk 4: The acceptance rejection method
####################################################

### Sampling from a beta distribution
f <- function(x) dbeta(x,shape1=3,shape2=3)
out <- optimize(f,interval=c(0,1),maximum=TRUE)
c <- out$objective
x <- seq(from=0,to=1,by=0.01)
plot(x,f(x),type='l')
set.seed(55555)
for (k in 1:1000){
  U1 <- runif(1)
  U2 <- runif(1,min=0,max=c)
  if (U2 <= f(U1))
    points(U1,U2,pch=20)
  else
    points(U1,U2,pch=1)
}


####################################################
#   Chunk 5: Sampling from the normal distribution
####################################################

### Interaction of the Box-Muller approach with
### a bad LCG
m <- 2048
a <- 1229
c <- 1
N <- m-2
seed <- 0
U <- LCG(a,c,m,seed,N)$U
U1 <- U[seq(from=1,by=2,to=N-1)]
U2 <- U[seq(from=2,by=2,to=N)]
X <- sqrt(-2*log(U1)) * cos(2*pi*U2)
Y <- sqrt(-2*log(U1)) * sin(2*pi*U2)
par(mfrow=c(2,1))
plot(X,Y,pch=20,cex=0.5)
X <- sqrt(-2*log(U2)) * cos(2*pi*U1)
Y <-sqrt(-2*log(U2)) * sin(2*pi*U1)
plot(X,Y,pch=20,cex=0.5)
par(mfrow=c(1,1))

### Sampling a multivariate normal by a homemade function
Sigma <- matrix(c(4,1,-2,1,3,1,-2,1,5),nrow=3)
mu <- c(8,6,10)
# check eigenvalues !
eigen(Sigma)$values
set.seed(55555)
Z <- MultiNormrnd(mu,Sigma,10000)
colMeans(Z)
cov(Z)

### Finding the square root matrix by diagonalization
eig <- eigen(Sigma)
V <- eig$vectors
V %*% diag(eig$values) %*% t(V)
sqrtSigma <- V %*% diag(sqrt(eig$values)) %*% t(V)
sqrtSigma %*% sqrtSigma

### Sampling a multivariate normal by a library function
library(MASS)
set.seed(55555)
Z <- mvrnorm(n=10000,mu=mu,Sigma=Sigma)
colMeans(Z)
cov(Z)

####################################################
#   Chunk 6: Dealing with copulas
####################################################

### Sampling from a multivariate t copula 
### with gamma and normal marginals
set.seed(55555)
mu <- c(0,0)
Sigma <- matrix(c(1,0.8,0.8,1), nrow=2)
nu <- 10
n <- 10000
Y <- rmvt(n,delta=mu,sigma=Sigma,df=nu)
U <- pt(Y,df=nu)
X <- cbind(qgamma(U[,1],shape=2,scale=1), qnorm(U[,2]))
plot(X,pch=20,cex=0.5)
par(mfrow=c(2,1))
hist(X[,1],nclass=100)
hist(X[,2],nclass=100)
par(mfrow=c(1,1))

### An illustration of the copula package 
tc <- tCopula(param=c(0.8), dim=2, dispstr="ex", df=10)
# the following statement is a possible alternative
# tc <- ellipCopula(param=c(0.8),family='t',dim=2,
#                   dispstr="ex",df=10)
set.seed(55555)
U <- rCopula(10000, tc)
X <- cbind(qgamma(U[,1],shape=2,scale=1), qnorm(U[,2]))
par(mfrow=c(1,2))
plot(X,pch=20,cex=0.5)
Uhat <- pobs(X)
plot(U,pch=20,cex=0.5)
par(mfrow=c(1,1))
fitted.copula <- fitCopula(copula=tCopula(dim = 2),data=Uhat)
fitted.copula
