rm(list=ls())
require(zoo)
require(ggplot2)
require(rugarch)
require(MASS)
require(moments)

source("LoadFundData.R")
# source("LoadWigData.R")

par(mfrow=c(2,1), cex = 0.7, bty="l")
plot(P, main="Level", xlab="", ylab="")
plot(r, main="returns", xlab="", ylab="")

#####################################
# VaR model settings                #
#####################################

T  <- length(r)                        # full sample, nobs
N  <- 1250                             # we shorten the sample to the last 5 years
r  <- tail(r,N)                        # log-returns
R  <- coredata(r)

p  <- 0.05                             # tolerance level 
H  <- 4                                # horizon

##########################################################
# Square root of time method for normal distribution     #
##########################################################

m   <- mean(R) # m=0
s   <- sd(R)
mH  <- m*H
sH  <- s*sqrt(H)

VaR1_N    <- qnorm(p)*s + m            # VaR for h=1
ES1_N     <- m - s*dnorm(qnorm(p))/p   # Expected shortfall for h=1

VaRH_N    <- mH + sH*qnorm(p)            # VaR for horizon H
ESH_N     <- mH - sH*dnorm(qnorm(p))/p   # Expected shortfall for horizon H


# density plot
x = seq(-4*sH,4*sH,0.01*sH)
df1 <- as.data.frame(cbind(x,dnorm(x,m,s))); names(df1) <- c("r","y"); df1$horyzont = "h=1" 
dfH <- as.data.frame(cbind(x,dnorm(x,mH,sH))); names(dfH) <- c("r","y"); dfH$horyzont = "h=H" 
df  = rbind(df1,dfH)

ggplot(data = df, aes(x = r, y=y, color=horyzont)) + 
  geom_line(size=1.2) + 
  scale_color_grey() + 
  theme_classic() +
  labs(title="Distribution of returns for longer horizons", y="", x="returns", caption=" ")+
  geom_vline(xintercept=c(VaRH_N,VaR1_N), linetype="dashed", size=1.2, color=c("gray","black"))  

###############################
# Cornish-Fisher expansion    #
###############################
require(moments)
m   = mean(R)
s   = sd(R)
S   = skewness(R)
K   = kurtosis(R)-3

PsiP    <- qnorm(p)
VaR1_CF  <- m + s*(PsiP + (PsiP^2-1)/6*S + (PsiP^3-3*PsiP)/24*(K) - (2*PsiP^3-5*PsiP)/36*(S^2) ) 

mH  <- m*H
sH  <- s*sqrt(H)
SH  <- S/sqrt(H)
KH  <- K/H

VaRH_CF        <- mH + sH*(PsiP + (PsiP^2-1)/6*SH + (PsiP^3-3*PsiP)/24*(KH) - (2*PsiP^3-5*PsiP)/36*(SH^2) ) 
# VaRHCF        <- m*H + s*sqrt(H)*(PsiP + (PsiP^2-1)/6*(S/sqrt(H)) + (PsiP^3-3*PsiP)/24*(K/H) - (2*PsiP^3-5*PsiP)/36*(S^2/H) ) 

round(-100*c(VaR1_N, VaR1_CF),2)
round(-100*c(VaRH_N, VaRH_CF),2)

##########################################################
# Mnte Carlo simulations for normal distribution         #
##########################################################

# number of simulations
M      <- 10000      

# A matrix with draws
draws  <- matrix(NA,H,M)
# each column for each simulation
for(i in 1:M){
  draws[,i] <- rnorm(H, mean=m, sd=s)
}

# faster method
# draws <- matrix(rnorm(M*H, mean=m, sd=s),nrow=H,ncol=M)

draws0      <- colSums(draws)               # cumulated growth 
draws0      <- sort(draws0)                 # sorted obs
M0          <- floor(M*p)                   # p-th quantile                   
VaRH_NMC    <- draws0[M0]                   # compare to: quantile(draws0,p)

# ES
ESH_NMC     <- mean(draws0[1:M0])           

# Coparing numerical and analytical methods
round(-100*c(VaRH_N, VaRH_NMC),2)
round(-100*c(ESH_N, ESH_NMC),2)

######################################
# Monte Carlo for t-Student          #
######################################
require(MASS)
M      <- 10000       
R0 <- (R-m)/s
v     <- 5                                  # degrees of freedom 
#dt    <- fitdistr(R0, "t", m = 0, start = list(s=sqrt((v-2)/v), df=v), lower=c(0.001,3))
#v     <- dt$estimate[[2]]                   

# MC simulations 
draws  <- matrix(NA,H,M)
for(i in 1:M){
  draws[,i] <- m + sqrt((v-2)/v)*rt(H, v)*s
}
draws1      <- colSums(draws) 
draws1      <- sort(draws1)                
M0          <- floor(M*p)                                   
VaRH_t      <- draws1[M0]                   
ESH_t       <- mean(draws1[1:M0])           

# plot 
df0 = as.data.frame(draws0); names(df0) = "draws"; df0$rozklad = "Normal"
df1 = as.data.frame(draws1); names(df1) = "draws"; df1$rozklad = "t-Student"
df  = rbind(df0,df1)

ggplot(df, aes(x=draws, color=rozklad)) +
  geom_density(size=1.2)+
  scale_color_grey() + 
  theme_classic() +
  labs(title="Distribution of returns for longer horizons", y="", x="returns", caption=" ") +
  geom_vline(xintercept=c(VaRH_t,VaRH_N), linetype="dashed", size=1.2, color=c("gray","black"))  

##############################################
# Monte Carlo for GARCH model                #
##############################################
require(rugarch)
M      <- 10000       

# Specyfication: GARCH(1,1) 
garch.spec <- ugarchspec(variance.model = list(model="sGARCH", garchOrder=c(1,1)),
                        mean.model = list(armaOrder=c(0,0), include.mean=TRUE),
                        distribution.model = "norm")
#estimation
garch.fit  <- ugarchfit(data=R, spec=garch.spec)

# MC simulations
garch.sim   <- ugarchsim(garch.fit, n.sim = H, n.start = 0, m.sim = M, startMethod = "sample")
draws       <- garch.sim@simulation$seriesSim

# VaR/ES
draws2      <- colSums(draws) 
draws2      <- sort(draws2)                
M0          <- floor(M*p)                                    
VaRH_GARCH  <- draws2[M0]                   
ESH_GARCH   <- mean(draws2[1:M0])           

# plot
df2 = as.data.frame(draws2); names(df2) = "draws"; df2$rozklad = "GARCH"
df  = rbind(df0,df2)

ggplot(df, aes(x=draws, color=rozklad)) +
  geom_density(size=1.2)+
  scale_color_grey() + 
  theme_classic() +
  labs(title="Distribution of returns for longer horizons", y="", x="returns", caption=" ") +
  geom_vline(xintercept=c(VaRH_N,VaRH_GARCH), linetype="dashed", size=1.2, color=c("gray","black"))  


##########################################################
# BOOTSTRAP - historical simulation                      #
##########################################################

require(MASS)
M      <- 10000       

# Bootstraped values 
draws  <- matrix(NA,H,M)
for(i in 1:M){
  draws[,i] <- sample(R, H, replace=TRUE)  # draws from historical returns
}
draws3      <- colSums(draws) 
draws3      <- sort(draws3)                 
M0          <- floor(M*p)                              
VaRH_boot   <- draws3[M0]                   
ESH_boot   <- mean(draws3[1:M0])          

# Wykres 
df3 = as.data.frame(draws3); names(df3) = "draws"; df3$rozklad = "Bootstrap"
df  = rbind(df0,df3)

ggplot(df, aes(x=draws, color=rozklad)) +
  geom_density(size=1.2)+
  scale_color_grey() + 
  theme_classic() +
  labs(title="Distribution of returns for longer horizons", y="", x="returns", caption=" ") +
  geom_vline(xintercept=c(VaRH_N,VaRH_boot), linetype="dashed", size=1.2, color=c("gray","black"))  


########################################################
# For non IID returns                                  #
#########################################################

# H-period returns (they are H-times less numerous...)

PH    <- tail(P,N)                # Portfolio values
PH    <- PH[seq(N,1,-H)]          # obs each H session
rH    <- diff(log(PH))            # H-period returns
length(rH)

mNH   <- mean(rH) 
sNH   <- sd(rH)

VaRH_NH    <- qnorm(p)*sNH + mNH            
ESH_NH     <- mNH - sNH*dnorm(qnorm(p))/p   

round(-100*c(VaRH_NH,VaRH_N),3)
round(-100*c(ESH_NH,ESH_N),3)

x = seq(-4*sH,4*sH,0.01*sH)
dfNH <- as.data.frame(cbind(x,dnorm(x,mNH,sNH))); names(dfNH) <- c("r","y"); dfNH$method = "H-period returns" 
dfH  <- as.data.frame(cbind(x,dnorm(x,mH,sH))); names(dfH) <- c("r","y"); dfH$method = "1-period returns" 
df  = rbind(dfNH,dfH)

ggplot(data = df, aes(x = r, y=y, color=method)) + 
  geom_line(size=1.2) + 
  scale_color_grey() + 
  theme_classic() +
  labs(title="VaR for longer horizons", y="", x="st. zwrotu", caption=" ")+
  geom_vline(xintercept=c(VaRH_NH,VaRH_N), linetype="dashed", size=1.2, color=c("gray","black"))