En esta oportunidad seguiremos practicando análisis estadÃstico de series de tiempo. Primero, debemos cargar las librerias.
## [[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"
El paquete pdfetch nos permite descargar series temporales acciones cotizadas en mercados internacionales. En el siguiente código, obtenemos los datos de del Neovasc.
Neovasc <- pdfetch_YAHOO("NVCN",from = as.Date("2019-01-01"),to = as.Date("2020-01-01"), interval = '1d') #DATOS DE Neovasc
Podemos ver un resumen de sus estadÃsticas con la función summary. A este tipo de conjunto de datos se les conoce como Open-high-low-close chart (OHLC) permite analizar los diferentes movimientos de un instrumento financiero.
summary(Neovasc)
## Index NVCN.open NVCN.high NVCN.low
## Min. :2019-01-02 Min. : 2.570 Min. : 2.710 Min. :2.340
## 1st Qu.:2019-04-02 1st Qu.: 3.380 1st Qu.: 3.460 1st Qu.:3.203
## Median :2019-07-02 Median : 4.295 Median : 4.385 Median :4.030
## Mean :2019-07-02 Mean : 4.494 Mean : 4.663 Mean :4.277
## 3rd Qu.:2019-10-01 3rd Qu.: 5.125 3rd Qu.: 5.325 3rd Qu.:4.900
## Max. :2019-12-31 Max. :10.000 Max. :11.000 Max. :7.700
## NVCN.close NVCN.adjclose NVCN.volume
## Min. :2.580 Min. :2.580 Min. : 4100
## 1st Qu.:3.297 1st Qu.:3.297 1st Qu.: 52000
## Median :4.185 Median :4.185 Median : 107650
## Mean :4.431 Mean :4.431 Mean : 1071259
## 3rd Qu.:5.100 3rd Qu.:5.100 3rd Qu.: 807575
## Max. :8.300 Max. :8.300 Max. :61529000
Ahora crearemos un objeto Time-Series, el cual indica al software que se tomaron muestras equidistantes en el tiempo (meses, dias, años).
tsNeovasc <- ts(Neovasc$`NVCN.close`,start = c(2019,1),frequency=356.25)
*(Discretos)
\[r_t=(P_t-P_{t-1})/P_{t-1}\]
d_Neovasc <- diff(tsNeovasc)/tsNeovasc[-length(tsNeovasc)]
mu <- mean(d_Neovasc)
s2 <- var(d_Neovasc)
s <- sd(d_Neovasc)
*Retornos (Continuo) \[R_t=\text{ln}(P_t/P_{t-1})\]
Se puede demostrar que:
\[R_t=\text{ln}(1+r_t)\] Podemos calcular los retornos por medio de la diferencia de logaritmos con la función diff:
l_Neovasc<-diff(log(tsNeovasc))
muC<-mean(l_Neovasc)
s2C<-var(l_Neovasc)
sC<-sd(l_Neovasc)
Como se puede observar tanto mi media tienden a parecer a acercarse a cero, ya que está en un entorno de -0.0006067751, y esa media ha sido obtenida de un Retorno continuo, que en este caso este último se eléboro con el precio de periodo entre el precio del periodo anterior. Mientras que la varianza, que es una medida de dispersión que representa la variabilidad de una serie de datos respecto a su media. Formalmente se calcula como la suma de los residuos al cuadrado divididos entre el total de observaciones. En este caso se obtuvo 0.003433411. En este caso de la dispersión estándar se obtuvo 0.05859531, en este caso se puede decir que la dispersión de los datos alrededor de la media es cercano de cero.
muC
## [1] -0.0006067751
s2C
## NVCN.close
## NVCN.close 0.003433411
sC
## [1] 0.05859531
Como se puede observar el siguiente grafico en donde se puede ver el eje Y que mostraran la densidad de la distribución, también tenemos el eje X en las que está representando la expresión de los retornos continuos que son menores a cero. Dentro de la gráfica podemos ver dos lÃneas de distribución de una normal que estan representados por la media y la varianza de color rojo. Mientras que en la lÃnea punteada de color negro es la distribución de los retornos continuos. En el gráfico podemos observar que la forma de histograma no se asemeja a una forma asimétrica, por lo que nos podria decir que no tendria una distribucion normal, mientras que se puede ver que la distribución de los retornos continuos se aleje de una distribución normal, y la distribucion normal es graficada de manera muy alta propia de la varianza.
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)
5. Realice tres test de normalidad 1) Jarque-Bera, 2) Kolmogorov, 3) Spahiro-Wilk, ¿Cuál es la hipótesis nula de cada uno de estos test? Interprete sus resultados.(3pts)
Podemos realizar siguientes test de normalidad:
H0: La distribución es normal, H1: La distribución no es normal
A un nivel de significancia de 0.05, podemos ver que el p-value es menor al nivel de significancia lo que quiere decir que se rechaza la hipótesis nula y se acepta la hipótesis que afirma la no normalidad de los datos.
H0: La distribución es normal, H1: La distribución no es normal
Comprobamos el nivel de significación, si es menor que 0.05 la distribución no es normal, por lo que quiere decir que la distribución no es normal.
H0 : curtosis y asimetria de una normal
En este caso Jarque Bera nos muestra una probabilidad muy pequeña, menor al 0.05, por lo que nos quiere decir que se rechaza la hipotesis nula lo que quiere decir que su asimetria y su curtosis tienen comportamientos no normales.
shapiro.test(l_Neovasc)
##
## Shapiro-Wilk normality test
##
## data: l_Neovasc
## W = 0.8797, p-value = 0.0000000000003437
ks.test(l_Neovasc, "pnorm", mean=mu, sd=s)
## Warning in ks.test(l_Neovasc, "pnorm", mean = mu, sd = s): ties should not be
## present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: l_Neovasc
## D = 0.14717, p-value = 0.00003793
## alternative hypothesis: two-sided
jarque.bera.test(l_Neovasc)
##
## Jarque Bera Test
##
## data: l_Neovasc
## X-squared = 815.71, df = 2, p-value < 0.00000000000000022
Analizando la serie temporal:
tsNeovasc %>% ggtsdisplay()
Analizando la serie temporal de la diferencia de logaritmos con 1 rezago:
log(tsNeovasc) %>% diff(lag=1) %>% ggtsdisplay()
6. Para un nivel de significancia de α1 = 0.01, α2 = 0.05, α3 = 0.1 Calcule el Value at Risk, con los siguientes m´etodos (Exprese sus resultados en terminos monetarios para para una cartera de 100 000 d´olares).(6pts) a) Value at Risk de con datos hist´oricos para cada α b) Value at Risk asumiendo una que los datos se distribuyen como una normal con media µ y varianza σ^2, para un nivel de confianza de α. c) Value at Risk - Montecarlo Para cada α realice las siguiente simulaci´on: Realice 10 000 simulaciones de un proceso browniano con media µ y varianza σ^2 (A partir de los retornos continuos de la pregunta 2) Calcule el VaR de cada simulaci´on y guarde su resultado en un vector de datos. Grafique la distribuci´on de los VaR simulados. Calcule el promedio, y su intervalo de confianza con los percentiles 0.025 y 0.975 de los datos simulados.
Se define como el montó mÃnimo que se perderÃa su un evento adverso ocurriera en un determinado horizonte de tiempo.
riesgo estadistico 1) Calculamos la media y la varianza de los retornos
En el primer cálculo del VAR llamado VAR1 que es utilizado para medir el riesgo financiero, en este caso de Neovasc. Este caso podemos observar que entre el periodo del primero de enero del 2019 y el primero de enero del 2020, la máxima perdida serÃa de $-12,746,617.00 dólares, dado un nivel de confianza de 0.01.
W<-100000
alpha1 <- 0.01
q1 <- quantile(x=l_Neovasc, alpha1)
mean(exp(q1)-1)
## [1] -0.1274662
VAR1 <- W*(exp(q1)-1)
VAR1_Moneda<-(VAR1*W)/100
VAR1_Moneda
## 1%
## -12746617
Para un nivel de confianza de 0.05, se obtuvo que la máxima perdida de al redor de $7,260,825.00 dólares dúrate un año.
W<-100000
alpha2 <- 0.05
q2 <- quantile(x=l_Neovasc, alpha2)
mean(exp(q2)-1)
## [1] -0.07260825
VAR2 <- W*(exp(q2)-1)
VAR2_Moneda<-(VAR2*W)/100
VAR2_Moneda
## 5%
## -7260825
En el VAR3 podemos observar la medición en función al valor del portafolio que es de 100000 dólares, en donde se puede obtener una pérdida de $5,315,615.00 a un nivel de confianza del 0.1.
W<-100000
alpha3 <- 0.1
q3 <- quantile(x=l_Neovasc, alpha3)
mean(exp(q3)-1)
## [1] -0.05315615
VAR3 <- W*(exp(q3)-1)
VAR3_Moneda<-(VAR3*W)/100
VAR3_Moneda
## 10%
## -5315615
Con la media y la varianza o su raiz cuadrada que vendrÃa a ser la desviación estandar, podemos obtener mucha información de esta variable aleatoria. Por ejemplo, podemos hacer un gráfico de la distribución de los datos, y ajustar una distribución normal \(X \thicksim N(\mu=\) 0.0011 \(,\sigma=\) 0.06 \()\). Donde \(\mu\) es la media, y \(\sigma\) es la desviación estandar.
Podemos, por ejemplo, calcular la siguiente pregunta:
En un lenguaje más estadÃstico, nos piden calcular el valor crÃtico de nuestra distribución \(X \thicksim N(\mu=-0.0006,\sigma=0.06)\) correspondiente a una probabilidad del 1%, 5% y 10%.
set.seed(100000)
VARp1<-qnorm(0.01,mean = muC,sd = sC,lower.tail = T);VARp1
## [1] -0.1369199
#En términos continuos
exp(VARp1)-1
## [1] -0.1279599
set.seed(100000)
VARp2<-qnorm(0.05,mean = muC,sd = sC,lower.tail = T);VARp2
## [1] -0.09698749
#En términos continuos
exp(VARp2)-1
## [1] -0.09243264
set.seed(100000)
VARp3<-qnorm(0.1,mean = muC,sd = sC,lower.tail = T);VARp3
## [1] -0.07569969
#En términos continuos
exp(VARp3)-1
## [1] -0.07290542
El siguiente código nos permite graficar la distribución de nuestros datos y compararla con una\(X \thicksim N(\mu=\) 0.0011 \(,\sigma=\) 0.06 \()\). De izquierda a derecha, el área bajo la curva roja hasta la linea roja \((x=\)-0.1369\()\) es 1%, que és la probabilidad que nuestros retornos caigan en un punto de esa área.
x<-seq(-0.2,0.2,by=0.001)
plot(density(l_Neovasc),col="blue",main = c("Distribucion de los retornos de IBM"))
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 \(10\ 000\) observaciones de una \(X \thicksim N(0.001,0.064)\), y obtener los retornos esperados en el percentil 1
set.seed(1000)
simulacion <- rnorm(10000,mean = muC,sd = sC)
VARsim<-quantile(simulacion,0.01);VARsim
## 1%
## -0.1361785
exp(VARsim)-1
## 1%
## -0.1273132
VAR.mc1 <- numeric()
for (i in 1:10000) {
changes <- rnorm(length(l_Neovasc),mean=1+muC,sd=sC)
sim.ts <- cumprod(c(as.numeric(tsNeovasc[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.mc1[i] <- sim.VAR
}
mean(VAR.mc1)
## [1] -0.1324012
sd(VAR.mc1)
## [1] 0.01228368
plot(density(VAR.mc1))
quantile(VAR.mc1,0.025)
## 2.5%
## -0.1581742
quantile(VAR.mc1,0.975)
## 97.5%
## -0.1100444
VAR.mc2 <- numeric()
for (i in 1:10000) {
changes <- rnorm(length(l_Neovasc),mean=1+muC,sd=sC)
sim.ts <- cumprod(c(as.numeric(tsNeovasc[1]),changes))
sim.R <- diff(log(sim.ts))
sim.q <- quantile(sim.R,0.05,na.rm = T)
sim.VAR <- exp(sim.q)-1
VAR.mc2[i] <- sim.VAR
}
mean(VAR.mc2)
## [1] -0.09579935
sd(VAR.mc2)
## [1] 0.007736851
plot(density(VAR.mc2))
quantile(VAR.mc2,0.025)
## 2.5%
## -0.1116297
quantile(VAR.mc2,0.975)
## 97.5%
## -0.08141372
VAR.mc3 <- numeric()
for (i in 1:10000) {
changes <- rnorm(length(l_Neovasc),mean=1+muC,sd=sC)
sim.ts <- cumprod(c(as.numeric(tsNeovasc[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.07497792
sd(VAR.mc3)
## [1] 0.006212608
plot(density(VAR.mc3))
quantile(VAR.mc3,0.025)
## 2.5%
## -0.08717469
quantile(VAR.mc3,0.975)
## 97.5%
## -0.06294227
#install.packages("rugarch")
require(rugarch)
## Loading required package: rugarch
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:purrr':
##
## reduce
## The following object is masked from 'package:stats':
##
## sigma
spec = ugarchspec()
fit = ugarchfit(data = Neovasc[,1], spec = spec)
fit
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 6.36108 0.511940 12.425 0.000000
## ar1 0.98605 0.011848 83.222 0.000000
## ma1 -0.17990 0.081625 -2.204 0.027527
## omega 0.00426 0.000137 31.133 0.000000
## alpha1 0.00000 0.000359 0.000 1.000000
## beta1 0.97966 0.005090 192.478 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 6.36108 0.193956 32.79642 0.00000
## ar1 0.98605 0.013887 71.00316 0.00000
## ma1 -0.17990 0.144232 -1.24728 0.21230
## omega 0.00426 0.004620 0.92203 0.35651
## alpha1 0.00000 0.000968 0.00000 1.00000
## beta1 0.97966 0.014613 67.04164 0.00000
##
## LogLikelihood : -171.2301
##
## Information Criteria
## ------------------------------------
##
## Akaike 1.4066
## Bayes 1.4906
## Shibata 1.4055
## Hannan-Quinn 1.4404
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.004954 0.9439
## Lag[2*(p+q)+(p+q)-1][5] 2.371163 0.8418
## Lag[4*(p+q)+(p+q)-1][9] 4.280436 0.6242
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1903 0.6627
## Lag[2*(p+q)+(p+q)-1][5] 0.2874 0.9850
## Lag[4*(p+q)+(p+q)-1][9] 0.3597 0.9995
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.01372 0.500 2.000 0.9068
## ARCH Lag[5] 0.04021 1.440 1.667 0.9964
## ARCH Lag[7] 0.08237 2.315 1.543 0.9996
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 7.4394
## Individual Statistics:
## mu 0.1721
## ar1 0.1569
## ma1 0.1157
## omega 0.1304
## alpha1 0.4879
## beta1 0.1299
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.7172 0.08719 *
## Negative Sign Bias 0.3761 0.70714
## Positive Sign Bias 0.9145 0.36135
## Joint Effect 6.6464 0.08406 *
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 190.7 0.000000000000000000000000000002388
## 2 30 209.2 0.000000000000000000000000000034209
## 3 40 214.7 0.000000000000000000000000039140218
## 4 50 238.5 0.000000000000000000000000010120781
##
##
## Elapsed time : 1.061063