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

#  Required packages
require(MASS)
require(lmtest)
require(forecast)

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

####################################################
#   Chunk 1: Confidence intervals for expected values 
####################################################
# Example 1
set.seed(55555)
X = rnorm(20,mean=10,sd=20)
t.test(X)
# Example 2
ci <- t.test(X,conf.level = 0.99)
ci$estimate
ci$conf.int

# Checking actual coverage 
MakeExpSample <- function(n,ev){rexp(n,rate=1/ev)}
MakeNormSample <- function(n,ev){rnorm(n,mean=ev,sd=ev)}
set.seed(55555)
nRuns <- 10000
sampleSize <- 20
alpha <- 0.05
ev <- 10
coverageExp <- CheckCoverage(ev, nRuns, SampleSize, 
                             MakeExpSample)
coverageExp 
coverageNorm <- CheckCoverage(ev, nRuns, SampleSize, 
                              MakeNormSample)
coverageNorm


####################################################
#   Chunk 2: Hypothesis testing
####################################################

### Testing the expected value
# two-tails
set.seed(55555)
X <- rnorm(10,10,20)
t.test(X,mu=20)
# one-tail
t.test(X,mu=20,alternative='less')
t.test(X,mu=20,alternative='greater')

### Correlation testing
set.seed(5)
Z <- rnorm(20)
X1 <- rnorm(20,10,5)
X2 <- rnorm(20,30,8)
cor.test(X1+5*Z,X2+5*Z)


####################################################
#   Chunk 3: Parameter estimation by MLE
####################################################

### Fitting a gamma distribution by MLE
set.seed(55555)
shape <- 2
rate <- 10
X <- rgamma(500, shape=shape, rate=rate)
# ---
library(stats4)
ll <- function(shape,rate) -sum(log(dgamma(X,shape,rate)))
fit1 <- mle(ll, start=list(shape=3,rate=10), nobs=length(X), 
            method = "L-BFGS-B", lower=c(0,0))
fit1@coef
confint(fit1)
# ---
fit2 <- fitdistr(X,"gamma")
fit2$estimate
fit2$sd


####################################################
#   Chunk 4: Checking distribution fit
####################################################

### Visual check by histograms and QQ-plots
set.seed(55555)
par(mfrow=c(1,2))
X=rnorm(200,10,40)
hist(X)
hist(X,nclass=50)
par(mfrow=c(1,1))
qqnorm(X)
qqline(X)

### Chi-square test
set.seed(55555)
shape <- 2
rate <- 10
X <- rgamma(500, shape, rate)
breakPts <- c(seq(0,0.8,0.1), 2)
probs <- pgamma(breakPts,shape,rate)
binProbs <- diff(probs)
f.obs <- hist(X, breakPts, plot = FALSE)$counts
chisq.test(x=f.obs,p=binProbs/sum(binProbs))
### Kolomogorov-Smirnov test
ks.test(X,pgamma, shape, rate)

### Testing normality 
set.seed(55555)
X=rnorm(20,10,20)
shapiro.test(X)
X=rexp(20,10)
shapiro.test(X)
X=rgamma(20,2,10)
shapiro.test(X)


####################################################
#   Chunk 5: Fitting linear regression models
####################################################

### Ordinary least squares, including testing residuals
set.seed(55555)
X1 <- 1:10
X2 <- rnorm(10,0,10)
eps <- rnorm(10,0,20)
Y <- 30 + 5*X1 + 3*X2 + eps
lmod1 <- lm(Y~X1+X2)
summary(lmod1)
confint(lmod1)
shapiro.test(lmod1$resid)
dwtest(lmod1)

### Polynomial regression
set.seed(55555)
X <- rnorm(10,0,10)
eps <- rnorm(10,0,20)
Y <- 300 + 5*X + 0.5*X^2 + eps
lmod2 <- lm(Y~poly(X,2,raw=TRUE))
summary(lmod2)

### Using interaction variables
X1 <- runif(40,min=0,max=50)
X2 <- runif(40,min=0,max=50)
Y1 <- 70 + 2*X1 + rnorm(40,0,10)
Y2 <- 30 + 5*X2 + rnorm(40,0,10)
Y <- c(Y1,Y2)
X <- c(X1,X2)
Z <- c(rep(1,40),rep(0,40))
lmod3 <- lm(Y~X*Z)
summary(lmod3)
plot(X1,Y1,pch=1,col="black",xlim=c(0,max(X)),
      ylim=c(0.8*min(Y),max(Y)),xlab="X",ylab="Y")
points(X2,Y2, pch=2,col="red")
abline(lmod3$coeff[1]+lmod3$coeff[3],
        lmod3$coeff[2]+lmod3$coeff[4],col="red")
abline(lmod3$coeff[1],lmod3$coeff[2],col="black")


####################################################
#   Chunk 6: Fitting time series models
####################################################

### Simulating and fitting an ARMA(2,2) model
set.seed(55555)
x1 <- arima.sim(n=1500, model=list(ar=c(0.85,-0.45), 
                                   ma=c(-0.23,0.25)), n.start=50, 
                rand.gen=function(n) 0.5*rt(n, df=4))
arima(x1,order=c(2,0,2))
auto.arima(x1)
par(mfrow=c(1,2))
ts.plot(x1)
ts.plot(x1[1:100])
par(mfrow=c(1,1))
x2 <- arima.sim(n=1500, list(ar=c(0.85,-0.45), 
                             ma=c(-0.23,0.25)), n.start=50, 
                rand.gen=function(n) 0.5*rt(n,df=4))
arima(x2,order=c(2,0,2))
auto.arima(x2)

####################################################
#   Chunk 7: Bayesian learning for coin flipping
####################################################

# BAYESIAN COIN FLIPPING
theta <- 0.2
# build a sample path of 1000 unfair coin flips, 
# fixing the initial four flips
set.seed(55555)
# the initial portion of the sample path is given
# the rest is random
history <- c(1, 1, 0, 0, ifelse(runif(996)<=theta,1,0))
# cumulate the number of heads at each time step
cumHistory <- c(0, cumsum(history))
# posterior likelihood, when h heads are drawn in n flips
post <- function(x,n,h) (x^h) * ((1-x)^(n-h))
# the nine time instants at which posterior is plotted
times <- c(0:4, 10, 20, 100, 1000)
# x-values to plot posterior densities
x <- seq(from=0, by=0.001, to=1)
# note the use of letters to select nine characters
# to be used later in printing labels for plots
labels <- letters[1:9]
par(mfrow=c(3,3))
for (i in 1:9){
  nFlips <- times[i]
  nHeads <- cumHistory[times[i]+1]
  y <- post(x, nFlips, nHeads)
  # sprintf is used to print a formatted string 
  lab <- sprintf('(%s) N=%d; H=%d', labels[i],nFlips,nHeads)
  plot(x,y/max(y),xlim=c(0,1),ylim=c(0,1.1),type='l',ylab='',
       xlab='',main=lab)
}
par(mfrow=c(1,1))
