#Bibliotecas
#---------------------------------------
library(pdfetch)
## Warning: package 'pdfetch' was built under R version 4.3.3
library(stats)
library(graphics)
library(grDevices)
library(utils)
library(datasets)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forcats)
library(stringr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(purrr)
library(tidyr)
library(forecast)
library(tibble)
library(ggplot2)

El paquete pdftech nos permite descargar series temporales acciones cotizadas en mercados internacionales.

#Datos
#------------------------------------------
Neovasc <- pdfetch_YAHOO("ALNEV.pA",from = as.Date("2019-01-01"),to = as.Date("2020-01-01"), interval = '1d')
head(Neovasc)
##            ALNEV.pA.open ALNEV.pA.high ALNEV.pA.low ALNEV.pA.close
## 2019-01-02      16240000      16760000     15680000       16720000
## 2019-01-03      17400000      17400000     16560000       16880000
## 2019-01-04      16640000      17160000     16560000       16640000
## 2019-01-07      16840000      17000000     16240000       16720000
## 2019-01-08      17400000      17720000     16280000       16600000
## 2019-01-09      16600000      17160000     16320000       16480000
##            ALNEV.pA.adjclose ALNEV.pA.volume
## 2019-01-02          16720000               0
## 2019-01-03          16880000               0
## 2019-01-04          16640000               0
## 2019-01-07          16720000               0
## 2019-01-08          16600000               0
## 2019-01-09          16480000               0
# Estadísticos descriptivos
#----------------------------------
summary(Neovasc)
##      Index            ALNEV.pA.open      ALNEV.pA.high       ALNEV.pA.low     
##  Min.   :2019-01-02   Min.   : 4120000   Min.   : 4400000   Min.   : 3920000  
##  1st Qu.:2019-04-01   1st Qu.: 7154000   1st Qu.: 7192000   1st Qu.: 6892000  
##  Median :2019-07-03   Median :10536000   Median :10688000   Median :10248000  
##  Mean   :2019-07-02   Mean   :11589250   Mean   :11867406   Mean   :11218344  
##  3rd Qu.:2019-10-01   3rd Qu.:16640000   3rd Qu.:16850000   3rd Qu.:16130000  
##  Max.   :2019-12-31   Max.   :27920000   Max.   :30560000   Max.   :25200000  
##  ALNEV.pA.close     ALNEV.pA.adjclose  ALNEV.pA.volume
##  Min.   : 4112000   Min.   : 4112000   Min.   :0      
##  1st Qu.: 7094000   1st Qu.: 7094000   1st Qu.:0      
##  Median :10376000   Median :10376000   Median :0      
##  Mean   :11502406   Mean   :11502406   Mean   :0      
##  3rd Qu.:16480000   3rd Qu.:16480000   3rd Qu.:0      
##  Max.   :26160000   Max.   :26160000   Max.   :0
# Rentabilidades
# --------------------------------------
ts <- ts(Neovasc$ALNEV.pA.close)
plot(ts)

# Retornos discretos
# -------------------------------
d_Neovasc <- diff(ts)/ts[-length(ts)]
mu <- mean(d_Neovasc)
s2 <- var(d_Neovasc)
s <- sd(d_Neovasc)
# Retoros continuos
#----------------------------------------
l_Neovasc <- diff(log(ts))
muC <- mean(l_Neovasc)
s2C <- var(l_Neovasc)
sC <- sd(l_Neovasc)
muC
## [1] -0.005027024
s2C
##                ALNEV.pA.close
## ALNEV.pA.close    0.001374129
sC
## [1] 0.03706925
CV = -sC/muC*100
CV
## [1] 737.3994

Histograma de los retornos continuos

x<-seq(-0.1,0.1,by=0.01)
hist(
     l_Neovasc,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_Neovasc),lwd=1.5,lty=2)
curve(dnorm(x,mean=muC,sd=s2C),lwd=2,lty=2,col="red",add = T)

Test de normalidad

#Contraste de normalidad Shapiro Wilks
#--------------------------------------
shapiro.test(l_Neovasc)
## 
##  Shapiro-Wilk normality test
## 
## data:  l_Neovasc
## W = 0.87143, p-value = 8.236e-14

El p valor 0.000 < 0.05. Rechazo H0. Los retornos no son normales.

# Contraste de normalidad Lilliefors
#-----------------------------------------
ks.test(l_Neovasc, "pnorm", mean = muC, sd = sC)
## Warning in ks.test.default(l_Neovasc, "pnorm", mean = muC, sd = sC): ties
## should not be present for the Kolmogorov-Smirnov test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  l_Neovasc
## D = 0.13794, p-value = 0.000122
## alternative hypothesis: two-sided

El p valor 0.0001 < 0.05. Rechazo H0. Los retornos no son normales.

#Test de Jarque Bera
#----------------------------
jarque.bera.test(l_Neovasc)
## 
##  Jarque Bera Test
## 
## data:  l_Neovasc
## X-squared = 608.41, df = 2, p-value < 2.2e-16

El p valor 0.000 < 0.05. Rechazo H0. Los retornos no son normales.

Análisis de la serie temporal

#Cotización
ts %>% ggtsdisplay()

En el primer gráfico, vemos que hay tendencia bajista. Y en los otros dos gráficos de autocorrelograma, las líneas de autocorrelograma traspasan los intervalos de confianza; por lo que vemos que la serie no es estacionaria.

Por tanto, hay que tomar una diferencia.

Analizando la serie temporal de la diferencia de logaritmos con 1 rezago (lag):

#Primer adiferencia del logaritmo ~ rentabilidad diaria 
#---------------------------------------------------------
log(ts) %>% diff(lag=1)  %>% ggtsdisplay()

Vemos que aquí, al no haber ningún tipo de tendencia en el primer gráfico y al estar las líneas autocorrelograma denteo de los intervalos de confianza,la serie sí es estacionaria.

Cuando el nivel es no estacionario y la primera diferencia es estacionaria y presenta un ruido blanco, tenemos un random walk (paseo aleatorio).

Valor en riesgo: con datos históricos

W <- 100000 # 100000 um
alpha1 <- 0.01
q1 <- quantile(x=l_Neovasc, alpha1)
mean(exp(q1)-1)
## [1] -0.09335893
VAR1 <- W*(exp(q1)-1)
VAR_um <- (VAR1*W/100)
VAR_um
##       1% 
## -9335893

La máxima pérdida esperada es 9.335.893 um

W <- 100000 # 100000 um
alpha2 <- 0.05
q2 <- quantile(x=l_Neovasc, alpha2)
mean(exp(q2)-1)
## [1] -0.05254355
VAR2 <- W*(exp(q2)-1)
VAR2_um <- (VAR2*W/100)
VAR2_um
##       5% 
## -5254355

La máxima pérdida esperada es 5.254.355 um

W <- 100000 # 100000 um
alpha3 <- 0.10
q3 <- quantile(x=l_Neovasc, alpha3)
mean(exp(q3)-1)
## [1] -0.0437085
VAR3 <- W*(exp(q3)-1)
VAR3_um <- (VAR3*W/100)
VAR3_um
##      10% 
## -4370850

La máxima pérdida esperada es 4.370.850 um

VAR paramétrico

set.seed(12345)
VARp1 <- qnorm(0.01,mean = muC, sd=sC, lower.tail = T)
VARp1
## [1] -0.091263
#En términos continuos
q4 <- exp(VARp1) - 1
q4
## [1] -0.08722238
VAR4 <- W*(exp(q4)-1)
VAR4_um <- (VAR4*W/100)
VAR4_um
## [1] -8352673

La máxima pérdida esperada es 8.352.673 um

set.seed(12345)
VARp2 <- qnorm(0.05,mean = muC, sd=sC, lower.tail = T)
VARp2
## [1] -0.06600051
#En términos continuos
q5 <- exp(VARp2) - 1
q5
## [1] -0.06386962
VAR5 <- W*(exp(q5)-1)
VAR5_um <- (VAR5*W/100)
VAR5_um
## [1] -6187269

La máxima pérdida esperada es 6.187.269 um

set.seed(12345)
VARp3 <- qnorm(0.10,mean = muC, sd=sC, lower.tail = T)
VARp3
## [1] -0.05253318
#En términos continuos
q6 <- exp(VARp3) - 1
q6
## [1] -0.05117716
VAR6 <- W*(exp(q6)-1)
VAR6_um <- (VAR6*W/100)
VAR6_um
## [1] -4988967

La máxima pérdida esperada es 4.988.967 um

x<-seq(-0.2,0.2,by=0.001)
plot(density(l_Neovasc),col="blue",main = c("Distribucion de los retornos de Neovasc"))
lines(x,dnorm(x,mu,s),col="red", lty = 1,lwd = 1)
lines(rep(VARp1,2),c(0,15),col="red",lwd = 1)
legend("topleft", legend = c("Datos reales", "Simulacion de una N(0.001,0.064)"),
       col = c("blue", "red"), lty = 1, cex = 0.8)

Podemos simular 1000000 observaciones de una X~N(0.001,0.0064), y obtener los retornos esperados en el percentil 1.

set.seed(12345)
simulacion1 <- rnorm(1000000,mean = muC, sd=sC)
VARsim1 <- quantile(simulacion1,0.01)
VARsim1
##          1% 
## -0.09133733
hist(simulacion1)

Al hacer tantas distribuciones, vemos que este histograma tiene forma de distribución normal.

q7 <- exp(VARsim1)-1
q7
##          1% 
## -0.08729022
VAR7 <- W*(exp(q7)-1)
VAR7_um <- (VAR7*W/100)
VAR7_um
##       1% 
## -8358891

La máxima pérdida esperada es 8.358.891 um

VAR por simulación de Montecarlo

set.seed(12345)
VAR.mc3 <- numeric()
for (i in 1:10000) {
  changes <- rnorm(length(l_Neovasc),mean=1+muC,sd=sC)
  sim.ts <- cumprod(c(as.numeric(ts[1]),changes))
  sim.R <- diff(log(sim.ts))
  sim.q <- quantile(sim.R,0.1,na.rm = T)
  sim.VAR <- exp(sim.q)-1
  VAR.mc3[i] <- sim.VAR
}
mean(VAR.mc3)
## [1] -0.05215722
plot(VAR.mc3)

hist(VAR.mc3)

Vemos que tiene pinta de campaniforme.

sd(VAR.mc3)
## [1] 0.003932378
plot(density(VAR.mc3))

q8 <- quantile(VAR.mc3,0.01)
q8
##          1% 
## -0.06136139
q81 <- exp(q8)-1
q81
##         1% 
## -0.0595167
q81 <- exp(q8)-1
q81
##         1% 
## -0.0595167
VAR8 <- W*(exp(q81)-1)
VAR8_um <- (VAR8*W/100)
VAR8_um
##       1% 
## -5778020

La máxima pérdida esperada es 5.778.020. um