#Instalando paquetes

## [[1]]
## [1] "pdfetch"   "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[2]]
## [1] "tseries"   "pdfetch"   "stats"     "graphics"  "grDevices" "utils"    
## [7] "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "forcats"   "stringr"   "dplyr"     "purrr"     "readr"     "tidyr"    
##  [7] "tibble"    "ggplot2"   "tidyverse" "tseries"   "pdfetch"   "stats"    
## [13] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[4]]
##  [1] "forecast"  "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
##  [7] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "tseries"   "pdfetch"  
## [13] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
## [19] "base"

#Importando datos

EXASdata <- pdfetch_YAHOO("^GSPC",from = as.Date("2019-01-01"),to = as.Date("2020-01-01"), interval = '1d')  #DATOS DE EXAS Exact Sciences Corporation

#Creando serie de tiempo

tsEXAS <- ts(EXASdata$`^GSPC.close`,start = c(2019,1),frequency=356.25)

##2: Calculando los retornos

1) Discretos

d_EXAS <- diff(tsEXAS)/tsEXAS[-length(tsEXAS)]
mu <- mean(d_EXAS)
s2 <- var(d_EXAS)
s <- sd(d_EXAS)

2) Continuos

l_EXAS<-diff(log(tsEXAS))
mu<-mean(l_EXAS)
s2<-var(l_EXAS)
s<-sd(l_EXAS)

##3: Análisis descriptivo ##4: Histograma de los retornos

x<-seq(-0.1,0.1,by=0.01)
hist(
     l_EXAS,prob=TRUE,ylim=c(0,80),xlim = c(-0.1,0.1),breaks = 50,col = "grey94",
     main = c("Histograma de los retornos"),
     xlab = expression(r==ln(P[t]/P[t-1])),
     ylab=c("Densidad"),
    )
lines(density(l_EXAS),lwd=1.5,lty=2)
curve(dnorm(x,mean=mu,sd=s),lwd=2,lty=2,col="red",add = T)

##6: Test de normalidad

1) Shapiro-Wilk

2) Kolmogorov

3) Jarque-Bera

shapiro.test(l_EXAS)
## 
##  Shapiro-Wilk normality test
## 
## data:  l_EXAS
## W = 0.94326, p-value = 0.00000002771
ks.test(l_EXAS, "pnorm", mean=mu, sd=s)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  l_EXAS
## D = 0.089311, p-value = 0.03648
## alternative hypothesis: two-sided
jarque.bera.test(l_EXAS)
## 
##  Jarque Bera Test
## 
## data:  l_EXAS
## X-squared = 125.76, df = 2, p-value < 0.00000000000000022

##6: Calcule el Value at Risk # a) Value at Risk de con datos historicos para cada

W<-100000

alpha <- 0.01
q1 <- quantile(x=l_EXAS, alpha)
mean(exp(q1)-1)
## [1] -0.02535171
VAR1 <- W*(exp(q1)-1)
W<-100000

alpha <- 0.05
q1 <- quantile(x=l_EXAS, alpha)
mean(exp(q1)-1)
## [1] -0.01208628
VAR1 <- W*(exp(q1)-1)
W<-100000

alpha <- 0.1
q1 <- quantile(x=l_EXAS, alpha)
mean(exp(q1)-1)
## [1] -0.00728274
VAR1 <- W*(exp(q1)-1)

b) Value at Risk asumiendo una que los datos se distribuyen como una normal con media

set.seed(1000)
VARp<-qnorm(0.01,mean = mu,sd = s,lower.tail = T);VARp
## [1] -0.01733375
#En términos discretos
exp(VARp)-1
## [1] -0.01718439

#c) Value at Risk - Montecarlo

set.seed(1000)
simulacion <- rnorm(10000,mean = mu,sd = s)
VARsim<-quantile(simulacion,0.01);VARsim
##          1% 
## -0.01723401
exp(VARsim)-1
##          1% 
## -0.01708635

#MONTECARLO

VAR.mc <- numeric()
for (i in 1:1000) {
  changes <- rnorm(length(l_EXAS),mean=1+mu,sd=s)
  sim.ts <- cumprod(c(as.numeric(tsEXAS[1]),changes))
  sim.R <- diff(log(sim.ts))
  sim.q <- quantile(sim.R,0.01,na.rm = T)
  sim.VAR <- exp(sim.q)-1
  VAR.mc[i] <- sim.VAR
}
mean(VAR.mc)
## [1] -0.01662756
sd(VAR.mc)
## [1] 0.00166096
plot(density(VAR.mc))

quantile(VAR.mc,0.225)
##       22.5% 
## -0.01775673
quantile(VAR.mc,0.975)
##       97.5% 
## -0.01365065