##################################################################
### Example 8.9: Parametric bootstrap
##################################################################

#source("~/Documents/University/Book/R/VaR.r")
#source("~/Documents/University/Book/R/ES.r")
#source("~/Documents/University/Book/R/VaRexactconfbounds.r")
source("~/Documents/University/Book/R/empquant.r")

numberofdays <- 500 # Selecting the size of the sample
data <- read.csv("Data/Nasdaq.csv", header = TRUE)
data <- data[3:(3+numberofdays),]
dates.nq <- as.Date(data[,1], "m/%d/%y")
index.nq <- rev(data[,7])
logret <- diff(log(index.nq)) # Calculating log returns
n <- numberofdays-1
u <- (1:n)/(n+1)

###
# Functions
###

lstfun<-function(params,sdatavec)
{
n<-length(sdatavec)
sum((params[1]+params[2]*qt((1:n)/(n+1),params[3])-sdatavec)^2)
}

mltfit<-function(datavec,initguess)
{
n<-length(datavec)
lmtll<-function(params){-sum(log(dt((datavec-params[1])/params[2],params[3])/params[2]))}
optim(initguess,lmtll)
}


lstlsfit<-function(datavec,initguess)
{
n<-length(datavec)
sdatavec<-sort(datavec,dec=F)
locscatfun<-function(params){sum((params[1]+params[2]*qt((1:n)/(n+1),params[3])-sdatavec)^2)}
optim(initguess,locscatfun)
}

lsnlsfit<-function(datavec,initguess)
{
n<-length(datavec)
sdatavec<-sort(datavec,dec=F)
locscanfun<-function(params){sum((params[1]+params[2]*qnorm((1:n)/(n+1))-sdatavec)^2)}
optim(initguess,locscanfun)
}


###
# Normal and Student t model
###

N <- 1000
NormMLEvec <- matrix(NA,N,2)
NormLSEvec <- matrix(NA,N,2)
tMLEvec <- matrix(NA,N,3)
tLSEvec <- matrix(NA,N,3)

for(i in 1:N){
	logret.b <- sample(logret,n,replace=TRUE)
	#pquants.b <- empquant(logret.b,u)
	NormMLEvec[i,] <- c(mean(logret.b), sd(logret.b))
	NormLSEvec[i,] <-  lsnlsfit(logret.b, c(mean(logret.b), sd(logret.b)))$par
	tMLEvec[i,] <-  mltfit(logret.b, c(mean(logret.b), sd(logret.b),3))$par
	tLSEvec[i,] <- lstlsfit(logret.b, c(mean(logret.b), sd(logret.b),3))$par
}	

###
# Value-at-Risk
###

NormMLE <- c(mean(logret), sd(logret))
NormLSE <-  lsnlsfit(logret, c(mean(logret), sd(logret)))$par
tMLE <-  mltfit(logret, c(mean(logret), sd(logret),3))$par
tLSE <- lstlsfit(logret, c(mean(logret), sd(logret),3))$par

p <- 0.05
ResNormMLE <- 100*(1-exp(NormMLE[1]+NormMLE[2]*qnorm(p))) - 100*(1-exp(NormMLEvec[,1]+NormMLEvec[,2]*qnorm(p)))
ResNormLSE <- 100*(1-exp(NormLSE[1]+NormLSE[2]*qnorm(p))) - 100*(1-exp(NormLSEvec[,1]+NormLSEvec[,2]*qnorm(p)))
RestMLE <- 100*(1-exp(tMLE[1]+tMLE[2]*qt(p,tMLE[3]))) - 100*(1-exp(tMLEvec[,1]+tMLEvec[,2]*qt(p,tMLEvec[,3])))
RestLSE <- 100*(1-exp(tLSE[1]+tLSE[2]*qt(p,tLSE[3]))) - 100*(1-exp(tLSEvec[,1]+tLSEvec[,2]*qt(p,tLSEvec[,3])))
q <- 0.95
INormMLE <- c()
INormLSE <- c()
ItMLE <- c()
ItLSE <- c()

INormMLE <- c(100*(1-exp(NormMLE[1]+NormMLE[2]*qnorm(p))) +empquant(ResNormMLE,(1-q)/2), 100*(1-exp(NormMLE[1]+NormMLE[2]*qnorm(p))) +empquant(ResNormMLE,(1+q)/2))
INormLSE <- c( 100*(1-exp(NormLSE[1]+NormLSE[2]*qnorm(p))) + empquant(ResNormLSE,(1-q)/2), 100*(1-exp(NormLSE[1]+NormLSE[2]*qnorm(p))) + empquant(ResNormLSE,(1+q)/2))
ItMLE <- c(100*(1-exp(tMLE[1]+tMLE[2]*qt(p,tMLE[3])))+ empquant(RestMLE,(1-q)/2), 100*(1-exp(tMLE[1]+tMLE[2]*qt(p,tMLE[3])))+ empquant(RestMLE,(1+q)/2))
ItLSE <- c(100*(1-exp(tLSE[1]+tLSE[2]*qt(p,tLSE[3]))) + empquant(RestLSE,(1-q)/2), 100*(1-exp(tLSE[1]+tLSE[2]*qt(p,tLSE[3]))) + empquant(RestLSE,(1+q)/2))
INormMLE
INormLSE
ItMLE
ItLSE


hist(ResNormMLE, breaks=50, col="black", xlab="",ylab="", yaxt="n",main="", xlim = c(-0.6,0.6)) 


hist(ResNormLSE, breaks=50,col="black", xlab="",ylab="", yaxt="n",main="", xlim = c(-0.6,0.6)) 

hist(RestMLE, breaks=50,col="black", xlab="",ylab="", yaxt="n",main="", xlim = c(-0.6,0.6)) 

hist(RestLSE, breaks=50,col="black", xlab="",ylab="", yaxt="n",main="", xlim = c(-0.6,0.6)) 

###
# Expected Shortfall
###
p <- 0.05

ESNormMLE <- 100*(1-pnorm(qnorm(p)-NormMLE[2])*exp(NormMLE[1]+NormMLE[2]^2/2)/p)
ESNormLSE <- 100*(1-pnorm(qnorm(p)-NormLSE[2])*exp(NormLSE[1]+NormLSE[2]^2/2)/p)
intMLE <- integrate((function(x){100*(1-exp(tMLE[1]+tMLE[2]*qt(x,tMLE[3])))}),0,p)
EStMLE <- intMLE$value/p
intLSE <- integrate((function(x){100*(1-exp(tLSE[1]+tLSE[2]*qt(x,tLSE[3])))}),0,p)
EStLSE <- intLSE$value/p

ESNormMLE
ESNormLSE
EStMLE
EStLSE

ESResNormMLE <- ESNormMLE - 100*(1-pnorm(qnorm(p)-NormMLEvec[,2])*exp(NormMLEvec[,1]+NormMLEvec[,2]^2/2)/p)
ESResNormLSE <- ESNormLSE - 100*(1-pnorm(qnorm(p)-NormLSEvec[,2])*exp(NormLSEvec[,1]+NormLSEvec[,2]^2/2)/p)


ESRestMLE <- c()
ESRestLSE <- c()
for(i in 1:N){
	intMLEtmp <- integrate((function(x){100*(1-exp(tMLEvec[i,1]+tMLEvec[i,2]*qt(x,tMLEvec[i,3])))}),0,p)
	ESRestMLE <- c(ESRestMLE, EStMLE-intMLEtmp$value/p)
	intLSEtmp <- integrate((function(x){100*(1-exp(tLSEvec[i,1]+tLSEvec[i,2]*qt(x,tLSEvec[i,3])))}),0,p)
	ESRestLSE <- c(ESRestLSE, EStLSE-intLSEtmp$value/p)
}

q <- 0.95
IESNormMLE <- c(ESNormMLE + empquant(ESResNormMLE,(1-q)/2), ESNormMLE + empquant(ESResNormMLE,(1+q)/2))
IESNormLSE <- c(ESNormLSE + empquant(ESResNormLSE,(1-q)/2), ESNormLSE + empquant(ESResNormLSE,(1+q)/2))
IEStMLE <- c(EStMLE + empquant(ESRestMLE,(1-q)/2), EStMLE + empquant(ESRestMLE,(1+q)/2))
IEStLSE <- c(EStLSE + empquant(ESRestLSE,(1-q)/2), EStLSE +empquant(ESRestLSE,(1+q)/2))
IESNormMLE
IESNormLSE
IEStMLE
IEStLSE


hist(ESResNormMLE, breaks=50, col="black", xlab="",ylab="", yaxt="n",main="",xlim=c(-1.2,1.2)) 

hist(ESResNormLSE, breaks=50,col="black", xlab="",ylab="", yaxt="n",main="",xlim=c(-1.2,1.2)) 

hist(ESRestMLE, breaks=50,col="black", xlab="",ylab="", yaxt="n",main="",xlim=c(-1.2,1.2)) 

hist(ESRestLSE, breaks=50,col="black", xlab="",ylab="", yaxt="n",main="", xlim=c(-1.2,1.2)) 

###
# Polynomial Normal (not part of the Example)
###


pquants <- empquant(logret.nq,u)
plot(qnorm(u), pquants, xlab="",ylab="")

x <- qnorm(u)
x2 <- x^2
x3 <- x^3

mm <- glm(pquants~1+x+x2+x3)
inter <- mm$coefficients[[1]]
c1 <- mm$coefficients[[2]]
c2 <- mm$coefficients[[3]]
c3 <- mm$coefficients[[4]]
lines(x,inter+c1*x+c2*x2+c3*x3)


p <- 0.05

polynomVaR <- 100*(1-exp(inter+c1*qnorm(p)+c2*qnorm(p)^2+c3*qnorm(p)^3))
polynomVaR

polynomint <- integrate((function(x){100*(1-exp(inter+c1*qnorm(x)+c2*qnorm(x)^2+c3*qnorm(x)^3))}),0,p)
polynomES <- polynomint$value/p
polynomES


# Bootstrapping
N <- 1000
polynormVaRvec <- matrix(NA,N,1)
pred.logret <- polynormVaRvec
inter.b <- matrix(NA,N,1)
c1.b <- matrix(NA,N,1)
c2.b <- matrix(NA,N,1)
c3.b <- matrix(NA,N,1)

for(i in 1:N){
	logret.b <- sample(logret.nq,n,replace=TRUE)
	pquants.b <- empquant(logret.b,u)
	mm.b <- glm(pquants.b~1+x+x2+x3)
	inter.b[i] <- mm.b$coefficients[[1]]
	c1.b[i] <- mm.b$coefficients[[2]]
	c2.b[i] <- mm.b$coefficients[[3]]
	c3.b[i] <- mm.b$coefficients[[4]]
	polynormVaRvec[i] <- 100*(1-exp(inter.b[i]+c1.b[i]*qnorm(p)+c2.b[i]*qnorm(p)^2+c3.b[i]*qnorm(p)^3))
	z <- rnorm(1)
	pred.logret[i] <- inter.b[i]+c1.b[i]*z+c2.b[i]*z^2+c3.b[i]*z^3
	}
	

#for(i in 1:N) lines(x,inter.b[i]+c1.b[i]*x+c2.b[i]*x2+c3.b[i]*x3)

ResVaR <- polynomVaR - polynormVaRvec

polynomESvec <- c()
for(i in 1:N){
	polynominttmp <- integrate((function(x){100*(1-exp(inter.b[i]+c1.b[i]*qnorm(x)+c2.b[i]*qnorm(x)^2+c3.b[i]*qnorm(x)^3))}),0,p)
	polynomESvec <- c(polynomESvec, polynominttmp$value/p)
}

ResES <- polynomES-polynomESvec


IVaR <- c()
IES <- c()

IVaR <- c()
IVaR <- c(polynomVaR + empquant(ResVaR,(1-q)/2), polynomVaR + empquant(ResVaR,(1+q)/2))
IES <- c()
IES <- c(polynomES + empquant(ResES,(1-q)/2), polynomES + empquant(ResES,(1+q)/2))
IVaR
IES

hist(ResVaR,breaks = 50, col="black", xlab="",ylab="", yaxt="n",main="", xlim=c(-1.2,1.2)) 

hist(ResES,breaks = 50, col="black", xlab="",ylab="", yaxt="n",main="", xlim=c(-1.2,1.2))