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

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

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

####################################################
#   Chunk 1: Functions to deal with some continuous
#            distributions (univariate)
####################################################

### Uniform distribution
a <- 5
b <- 15
# PDF and CDF
dunif(8,min=a,max=b)
punif(8,min=a,max=b)
dunif(3,min=a,max=b)
punif(3,min=a,max=b)
dunif(18,min=a,max=b)
punif(18,min=a,max=b)
# CDF and quantiles are mutually inverse
punif(9,min=a,max=b)
qunif(0.4,min=a,max=b)
# random variates
runif(5,min=a,max=b)

### Binomial and geometric distributions
# plot PMF of binomial
par(mfrow=c(2,1))
plot(dbinom(0:30,size=30,p=0.2),type='h',lwd=5)
plot(dbinom(0:30,size=30,p=0.4),type='h',lwd=5)
par(mfrow=c(1,1))
# dealing with factorials and binomial coefficients
dbinom(10,size=30,p=0.4)
factorial(30)/factorial(20)/factorial(10)*0.4^10*0.6^20
choose(30,10)*0.4^10*0.6^20
# the geometric distribution
set.seed(55555)
rgeom(10,prob=0.5)

### Exponential and Poisson distributions
# Empirical check of the meoryless property
set.seed(55555)
CheckMemory(1,1,0.5,1000000)
# Check sample size using CDF 
1-pexp(1,rate=1)
1-pexp(0.5,rate=1)
# Plot PDF of Erlang variables
lambda <- 1
x <- seq(0,10,0.01)
nValues <- c(1,2,5,10)
data <- matrix(0,nrow=length(nValues),ncol=length(x))
for (i in 1:length(nValues))
   data[i,] <- dgamma(x, shape=nValues[i], rate = lambda)
plot(x,data[1,],ylim=c(0,max(data)),ylab="",type='l')
for (i in 2:length(nValues))
   lines(x,data[i,])
# Plot PMF of Poisson variable
plot(dpois(0:30,lambda=5),type='h',lwd=5)

### Normal and related distributions (chi2, Student, lognormal)
# Plot PDF of normal
x <- seq(-6.5,6.5,by=0.01)
plot(x,dnorm(x,mean=0,sd=1),type="l",lty="dashed")
lines(x,dnorm(x,mean=0,sd=2),type="l",lty="dotdash")
# CDF and quantiles of normal
pnorm(10,mean=5,sd=3)
pnorm((10-5)/3)
qnorm(0.95,mean=5,sd=3)
z=qnorm(0.95);z
5+3*z
# The 3 sigmas rule
pnorm(1) - pnorm(-1) 
pnorm(2) - pnorm(-2) 
pnorm(3) - pnorm(-3) 
# The lognormal distribution
set.seed(55555)
sigma <- 2
X <- rlnorm(1000000, meanlog=-(sigma^2)/2, sd=sigma)
mean(X)
# The chi-square distribution
x <- seq(0,20,by=0.01)
plot(x,dchisq(x,df=4),type="l",lty="dashed")
lines(x,dchisq(x,df=8),type="l",lty="dotdash")
set.seed(55555)
rchisq(3, df=4.5, ncp=3.6)
# The Stident's t distribution
x <- seq(-4.5,4.5,by=0.01)
plot(x,dnorm(x),type="l",lty="solid")
lines(x,dt(x,df=1),type="l",lty="dotdash")
lines(x,dt(x,df=5),type="l",lty="dashed")

# The beta distribution
x <- seq(0.001,0.999,by=0.001)
par(mfrow = c(3,3))
plot(x,dbeta(x,2,1),type='l')
plot(x,dbeta(x,1,1),type='l')
plot(x,dbeta(x,1,2),type='l')
plot(x,dbeta(x,1,0.8),type='l')
plot(x,dbeta(x,0.8,0.8),type='l')
plot(x,dbeta(x,0.8,1),type='l')
plot(x,dbeta(x,8,3),type='l')
plot(x,dbeta(x,3,3),type='l')
plot(x,dbeta(x,3,8),type='l')
par(mfrow = c(1,1))

# The gamma distribution
dgamma(2, shape=3, rate=2)
dgamma(2, shape=3, scale=1/2)
x <- seq(0,7,by=0.01)
plot(x,dgamma(x,shape=.5,rate=1),type='l',
        ylim=c(0,1.2),ylab='')
lines(x,dgamma(x,shape=1,rate=1))
lines(x,dgamma(x,shape=2,rate=1))
lines(x,dgamma(x,shape=3,rate=1))

# Empirical distributions and kernel densities
set.seed(55555)
data <- c(rnorm(40,10,10), rnorm(30,75,20))
hist(data,nclass=15)
windows()
par(mfrow=c(1,2))
kd1 <- density(data, bw=2)
plot(kd1,main="")
kd2 <- density(data, bw=10)
plot(kd2, main="")
class(kd1)
par(mfrow=c(1,1))


####################################################
#   Chunk 2: Multivariate distributions 
#            Covariance and correlation
####################################################

# Marginals are not enough
set.seed(55555)
len <- 50
times <- 0:len
eps <- rnorm(10)
lims <- c(len*min(eps),len*max(eps))
plot(x=times,y=eps[1]*times,type='l',ylim=lims)
points(x=times,y=eps[1]*times)
for (k in 2:10){
  lines(x=times,y=eps[k]*times)
  points(x=times,y=eps[k]*times)
}
windows()
eps <- rnorm(300)
path <- (0:300)*c(0,eps)
plot(x=0:300,y=path,type='l')

# Generating correlated normals
set.seed(55555)
rho <- 0.8
Z1 <- rnorm(100000)
Z2 <- rnorm(100000)
X1 <- Z1
X2 <- rho*Z1 + sqrt(1-rho^2)*Z2
cor(Z1,Z2)
cor(X1,X2)
cov(X1,X2)
# Using Cov and Cor
X <- as.matrix(cbind(2*X1,-3*X2))
cov(X)
cor(X)
cov2cor(cov(X))

# The limits of Pearson's correlation
set.seed(55555)
X=runif(100000)
cor(X,X)
cor(X,X^2)
cor(X,X^4)

# Cholesky factors and the square root matrix
Sigma <- matrix(c(10,2,2,8),nrow=2,ncol=2)
A <- chol(Sigma)
t(A) %*% A
sqrtSigma <- sqrtm(Sigma)
sqrtSigma
sqrtSigma %*% sqrtSigma

# Multivariate normal and t
mu <- c(20,10)
sigma <- matrix(c(25,35,35,100), nrow=2)
nu <- 2
pts <- seq(-20,60,0.1)
n <- length(pts)
X <- rep(pts,n)
Y <- rep(pts,rep(n,n))
grid <- cbind(X,Y)
densNorm <- dmvnorm(grid, mu,sigma,log=TRUE)
Z1 <- array(densNorm, dim=c(n,n))
contour(pts,pts,Z1,levels=seq(-10,-4))
windows()
densT <- dmvt(grid,mu,sigma,df=nu,log=TRUE)
Z2 <- array(densT, dim=c(n,n))
contour(pts,pts,Z2,levels=seq(-10,-4))

# Trying different correlation coefficients
set.seed(55555)
tryCorrelation()


####################################################
#   Chunk 3: Time series models
####################################################

### Dealing with seasonality 
# Detecting seasonality with ACF
set.seed(55555)
t <- 1:40
factorS <- c(15, 30, -10, -35)
Y <- 100+10*t+factorS[t%%4+1]+rnorm(length(t),0,20)
par(mfrow=c(1,2))
plot(t,Y,type='l')
points(t,Y,pch=15)
acf(Y)
# Detrending a seasonal time series
mod <- lm(Y~t)
dY <- Y-mod$fitted
plot(t,dY,type='l')
points(t,dY,pch=15)
acf(dY)
par(mfrow=c(1,1))

### Moving-average processes
# Example 1
set.seed(55555)
t <- 1:100
epsVet <- rnorm(101)
Y <- 40+epsVet[2:101]+0.8*epsVet[1:100]
par(mfrow=c(1,2))
plot(t,Y,type='l')
points(t,Y,pch=15)
acf(Y)
# Example 2
set.seed(55555)
t <- 1:100
epsVet <- rnorm(102)
Y<-40+epsVet[3:102]+0.8*epsVet[2:101]+0.6*epsVet[1:100]
par(mfrow=c(1,2))
plot(t,Y,type='l')
points(t,Y,pch=15)
acf(Y)
# Example 3
set.seed(55555)
t <- 1:100
epsVet <- rnorm(101)
Y <- 40+epsVet[2:101]+0.8*epsVet[1:100]
par(mfrow=c(1,2))
plot(t,Y,type='l')
points(t,Y,pch=15)
acf(Y)
# Using Partial ACF
expr1 = expression(Y[t] == 8 + 0.8 %*% Y[t-1] + epsilon[t])
pacf(Y, main = expr1)

### ARMA and ARIMA processes}
# Example 1
set.seed(5)
epsVet <- c(0, rnorm(100)) # first value is dummy
Y <- c(4,numeric(100))
for (t in 2:101)
  Y[t] <- Y[t-1] + 0.05*epsVet[t]
dY <- diff(Y)
par(mfrow=c(2,2))
plot(1:100,Y[-1],type='l')
points(1:100,Y[-1],pch=15)
acf(Y)
plot(1:100,dY,type='l')
points(1:100,dY,pch=15)
acf(dY)
par(mfrow=c(1,1))
# Example 2
tPhi1 <- t(matrix(c(0.7,0.2,0.2,0.7),2,2))
c.vec <- c(-0.7, 1.3)
mu.vec <- numeric(2)
cov.mat <- matrix(c(1,0.5,0.5,1),2,2)
set.seed(55555)
T <- 100
Ypath <- matrix(0, T, 2)
Ypath[1,] <- c(1,1) # initial condition
epsPath <- mvrnorm(n=T, mu.vec, cov.mat)
for (t in 2:T)
  Ypath[t,] <- c.vec+Ypath[t-1,]%*%tPhi1+epsPath[t,]
plot.ts(Ypath)

### Modeling stochastic volatility
# Sample path of a GARCH(1,1) model
set.seed(5555)
mu <- 0.0113
alpha0 <- 0.00092 
alpha1 <- 0.086
beta1 <- 0.853 
eps <- rnorm(550)
r <- rep(0, 550)
sigma2 <- rep(0, 550) 
u <- rep(0, 550)
for (i in 2:550) {
  sigma2[i] <- alpha0 + alpha1*u[i-1]^2 + beta1*sigma2[i-1] 
  u[i] <- eps[i] * sqrt(sigma2[i])
  r[i] <- mu + u[i] 
}
par(mfrow=c(2,1))
plot(r[51:550], type='l', ylab="return")
plot(sqrt(sigma2[51:550]), type='l', ylab="volatility")
par(mfrow=c(1,1))


####################################################
#   Chunk 4: Continuous-time processes
####################################################

### Sample paths of a standard Wiener process
set.seed(5)
T <- 1; numSteps <- 1000; numRepl <- 1
Wpath <- simWiener(T,numSteps,numRepl)
plot(0:numSteps,Wpath,type='l')

### Sample paths of a GBM
S0 <- 50; mu <- 0.1; sigma <- 0.3; T <- 1
numSteps <- 1000; numRepl <- 1
set.seed(5)
path <- simGBM(S0,mu,sigma,T,numSteps,numRepl)
plot(0:numSteps, path, type='l')


####################################################
#   Chunk 5: Principal component analysis
####################################################

### A simple numerical example
rho <- 0.85
sigma <- matrix(c(1, rho, rho, 1),2,2)
mu <- c(0, 0)
set.seed(55555)
NN <- 2000
X <- mvrnorm(NN, mu, sigma)
PCA <- prcomp(X, center=TRUE, scale=FALSE, retx=TRUE)
summary(PCA)
print(PCA)
# Checking the details of PCA "manually"
m <- colMeans(X);m
Xc <- X - (matrix(1,NN,1) %*% m)
S <- cov(X);S
vv <- eigen(S)
eigval <- vv$values;eigval
sqrt(eigval)
PCA$sdev
eigvet <- vv$vectors;eigvet
A <- PCA$rotation;A
z <- PCA$x
cov(z)
head(z)
head(Xc %*% A)
# Plotting principal components
par(mfrow=c(1,2))
plot(X[,1], X[,2])
plot(z[,1], z[,2])
par(mfrow=c(1,1))
boxplot(z[,1], z[,2], names=c("Z1","Z2"))
par(mfrow=c(1,1))

####################################################
#   Chunk 6: Risk-neutral derivative pricing
####################################################

### Plotting call option prices 
seqS0 <- 30:70; K <- 50; r <- 0.08; sigma <- 0.4
seqT <- seq(from=2, to=0, by=-0.25)
M <- length(seqT)
prices <- matrix(0,nrow=length(seqT),ncol=length(seqS0))
for (t in 1:length(seqT))
  prices[t,] <- EuVanillaCall(seqS0,K,seqT[t],r,sigma)
plot(seqS0, prices[1,], type='l',ylim=c(0,max(prices[1,])))
for (t in (2:(M-1))) 
  lines(seqS0, prices[t,])
lines(seqS0, pmax(0,seqS0-K),lwd=2)


