Modelo Garch
Equipo
2022-07-07
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(astsa)
library(forecast)##
## Attaching package: 'forecast'## The following object is masked from 'package:astsa':
##
## gaslibrary(quantmod)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTRlibrary(ggplot2)
library(tidyr)
library(MASS)
library(nortest)
library(moments)
library(PerformanceAnalytics)##
## Attaching package: 'PerformanceAnalytics'## The following objects are masked from 'package:moments':
##
## kurtosis, skewness## The following object is masked from 'package:graphics':
##
## legendlibrary(rugarch)## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:stats':
##
## sigmalibrary(ggtext)
library(GeneralizedHyperbolic)##
## Attaching package: 'GeneralizedHyperbolic'## The following object is masked from 'package:rugarch':
##
## qnig——– Modelo Garch (1,1)———–
Objetivo: Aplicar el modelo GARCH(1,1) para pronosticar la volatilidad de un activo financiero. Actividad: Cada equipo realizará el pronóstico de la volatilidad de 2 series de tiempo para un horizonte de 10 días en el futuro. Entregable: Se debe entregar el archivo .R que contenga el código utilizado para realizar el pronóstico.
1. Gráfica de la serie en niveles y dar una descripción.
#Obtener información en Yahoo
Televisa<-getSymbols("TV", from="2017-07-01", src = "yahoo",
auto.assign = F)[,6]
#Gráfica de serie de datos
ggplot(Televisa, aes(x=1:nrow(Televisa), y=TV.Adjusted))+
geom_line(col="darkorchid3",
size=1)+
labs(title="Serie de precios Televisa",
y="Precio",
x="Observaciones",
caption = "Fuente Yahoo Finance") +
theme(plot.title = element_text(color = "orange", size = 18, hjust = 0.5))
# 2. Obtener los rendimientos logarítmicos de la serie # 3. Graficar la
series de rendimientos logarítmicos y dar una descripción.
#En esta linea de código se estiman las diferencias logarítmicas aplicando de manera anidada las funciones diff() y log()
tv.L1.return <- diff(log(Televisa))
#se almacenan las diferencias logarítminas en un data frame para poder graficarlas mediante ggplot2()
tv.L1.return.df <- as.data.frame(tv.L1.return)
ggplot(tv.L1.return.df, aes(x=1:nrow(tv.L1.return.df),y=tv.L1.return))+
geom_line(col="darkorchid3",
size=1)+
labs(title="Retornos Logarítmicos Televisa",
y="Rendimiento",
x="Observaciones",
caption = "Fuente Yahoo Finance") +
theme(plot.title = element_text(color = "orange", size = 18, hjust = 0.5))## Don't know how to automatically pick scale for object of type xts/zoo. Defaulting to continuous.## Warning: Removed 1 row(s) containing missing values (geom_path).
qqnorm(tv.L1.return$TV.Adjusted, pch = 1, col="violet", frame = FALSE)
qqline(tv.L1.return$TV.Adjusted, col = "steelblue", lwd = 2)
Podemos observar que no se comporta como una distribución normal y lo podemos ver más claro con la segunda gráfica. También podemos hacer más pruebas de normalidad para corroborarlo.
4. Aplicar la prueba de normalidad Jarque-Bera a la serie de rendimientos logarítmicos y dar una descripción del resultado.
#Se estima la prueba de normalidad Jarque-Bera
tv.L1.return<-na.omit(tv.L1.return)
tv.L1.return.df<-na.omit(tv.L1.return.df)
jarque.bera.test(tv.L1.return)##
## Jarque Bera Test
##
## data: tv.L1.return
## X-squared = 2872.8, df = 2, p-value < 2.2e-16jarque.bera.test(tv.L1.return.df$TV.Adjusted)##
## Jarque Bera Test
##
## data: tv.L1.return.df$TV.Adjusted
## X-squared = 2872.8, df = 2, p-value < 2.2e-16Como p-value < .05 , entonces la serie no sigue una distribución normal.
5. Estimar la kurtosis de la serie de rendimientos logarítmicos y dar una interpretación del resultado.
kurtosis(tv.L1.return)## [1] 7.397179kurtosis(tv.L1.return.df)## [1] 7.397179Debido a que la kurtosis es >3, la distribución muestral es leptocúrticas. La asimetría es casi nula y es no-normal, kurtosis mayor a 3 y varianza diferente de 0, los rendimientos son negativos. También sugiere presencia de colas pesadas.
6. Estimar la volatilidad diaria y anual de la serie de retornos.
Volatilidad diaria:
sd(tv.L1.return)## [1] 0.02876257Volatilidad anualizada: Se considera 252 días que es el promedio de días de traiding que tienen un año. con lo que se obtiene la desviación estándar que se puede esperar de los retornos anualizados
sqrt(252)*sd(tv.L1.return)## [1] 0.4565916La volatilidad diaria es de 3%, mientras que la anualizada es del 45%.
#7. Estimar mediante la técnica “rolling window” la volatilidad mensual sobre escala anual.
chart.RollingPerformance(R=tv.L1.return,
width = 22,
FUN = "sd.annualized",
scale = 252,
main = "Rolling One Month Volatility")
Volatilidad mensual
8. Definir la especificación del modelo GARCH(1,1)
** ——– Modelo GARCH(1,1)———- **
#definición de parametros
alpha <- 0.1
beta <- 0.8
omega <- var(tv.L1.return)*(1-alpha-beta)
#Se obtienen los errores del modelo
e <- tv.L1.return-mean(tv.L1.return)
e2 <- e^2
#se obtiene la varianza para cada observación
nobs <- length(tv.L1.return)
predvar <- rep(NA,nobs)
#se inicializa el proceso
predvar[1] <- var(tv.L1.return)
#se inicializa el loop para estimar las varianzas
for(t in 2:nobs){
predvar[t] <- omega+alpha*e2[t-1]+beta*predvar[t-1]
}
#La volatilidad es la raíz cuadrada de la varianza
predvol <- sqrt(predvar)
predvol <- xts(predvol, order.by = time(tv.L1.return))11. Obtener la volatilidad incondicional
#se comparan los resultados con la volatilidad incondicional
uncvol <- sqrt(omega/(1-alpha-beta))
uncvol <- xts(rep(uncvol, nobs), order.by = time(tv.L1.return))#se genera la gráfica
plot(predvol)
lines(uncvol, col = "red", lwd =2)
Podemos observar que la volatilidad del Televisa tiene reversión a la media, se observan periodos por encima de la media y periodos por debajo de la media.
9. Realizar el ajuste de la serie de rendimientos al modelo GARCH(1,1)
#Errores de predicción
# Compute the mean daily return
m <- mean(tv.L1.return)
# Define the series of prediction errors
e <- tv.L1.return - m
# Plot the absolute value of the prediction errors
par(mfrow = c(2,1),mar = c(3, 2, 2, 2))
plot(abs(e))
# Plot the acf of the absolute prediction errors
acf(abs(e))
#Implementación en R utilizando la librería “rugrach”
#Especificación del modelo que queremos utilizar utilizando la función “ugarchspec()”, para este modelo vamos a especificar una media constante e igual a cero, para la varianza se utiliza un modelo GARCH estándar
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "sGARCH"),
distribution.model = "norm")Ajuste del modelo:
# Se estima el modelo GARCH(1,1) o sGARCH utilizando la función “ugarchfit()” utilizando los argumentos definidos en la especificación del modelo.
garchfit <- ugarchfit(data = tv.L1.return, spec = garchspec)#Pronóstico de la volatilidad de los rendimientos futuros, se utiliza la función “ugarchforecast()” utilizando el objeto que contiene el ajuste del modelo.
garchforecast <- ugarchforecast(fitORspec = garchfit, n.ahead = 5)Resultados ajuste del modelo:
garchcoef <- coef(garchfit)
garchunvar <- uncvariance(garchfit)
garchmean <- fitted(garchfit)
garchcoef## mu omega alpha1 beta1
## -1.382547e-03 8.598256e-06 7.080587e-02 9.215096e-0112. Plantear la ecuación del modelo con los parámetros obtenidos
- Ecuación * 2=8.55x10-6+.008e2+.09^2
sqrt(garchunvar)## [1] 0.0334500313. Generar la gráfica de volatilidades estimadas al ajustar el modelo
garchvol <- sigma(garchfit)
plot(garchvol)
lines(uncvol, col = "red", lwd =2)
tail(garchvol,1)## [,1]
## 2022-07-06 0.0282150914. Realizar el pronóstico del modelo
ug_forecast = ugarchforecast(garchfit, n.ahead = 30)
ug_forecast##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 30
## Roll Steps: 0
## Out of Sample: 0
##
## 0-roll forecast [T0=2022-07-06]:
## Series Sigma
## T+1 -0.001383 0.02734
## T+2 -0.001383 0.02739
## T+3 -0.001383 0.02744
## T+4 -0.001383 0.02749
## T+5 -0.001383 0.02754
## T+6 -0.001383 0.02759
## T+7 -0.001383 0.02764
## T+8 -0.001383 0.02769
## T+9 -0.001383 0.02774
## T+10 -0.001383 0.02779
## T+11 -0.001383 0.02784
## T+12 -0.001383 0.02788
## T+13 -0.001383 0.02793
## T+14 -0.001383 0.02798
## T+15 -0.001383 0.02802
## T+16 -0.001383 0.02807
## T+17 -0.001383 0.02812
## T+18 -0.001383 0.02816
## T+19 -0.001383 0.02820
## T+20 -0.001383 0.02825
## T+21 -0.001383 0.02829
## T+22 -0.001383 0.02834
## T+23 -0.001383 0.02838
## T+24 -0.001383 0.02842
## T+25 -0.001383 0.02846
## T+26 -0.001383 0.02850
## T+27 -0.001383 0.02855
## T+28 -0.001383 0.02859
## T+29 -0.001383 0.02863
## T+30 -0.001383 0.0286715. Obtener las volatilidades pronosticadas.
sigma(garchforecast)## 2022-07-06
## T+1 0.02733729
## T+2 0.02738946
## T+3 0.02744114
## T+4 0.02749232
## T+5 0.02754302ugarch2 = ugarchspec(mean.model = list(armaOrder = c(2,1)))
ugarch2##
## *---------------------------------*
## * GARCH Model Spec *
## *---------------------------------*
##
## Conditional Variance Dynamics
## ------------------------------------
## GARCH Model : sGARCH(1,1)
## Variance Targeting : FALSE
##
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model : ARFIMA(2,0,1)
## Include Mean : TRUE
## GARCH-in-Mean : FALSE
##
## Conditional Distribution
## ------------------------------------
## Distribution : norm
## Includes Skew : FALSE
## Includes Shape : FALSE
## Includes Lambda : FALSE10. Obtener la tabla de coeficientes haciendo el analizar de los t-values de cada parámetro.
ugfit = ugarchfit(spec = ugarch2, data = tv.L1.return)
ugfit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(2,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu -0.001401 0.000703 -1.9931 0.046247
## ar1 -0.436223 0.287209 -1.5188 0.128804
## ar2 0.072669 0.030074 2.4164 0.015676
## ma1 0.461919 0.286836 1.6104 0.107312
## omega 0.000008 0.000007 1.2112 0.225813
## alpha1 0.069502 0.013865 5.0127 0.000001
## beta1 0.922740 0.010630 86.8016 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu -0.001401 0.001145 -1.22348 0.221148
## ar1 -0.436223 0.232446 -1.87667 0.060564
## ar2 0.072669 0.026055 2.78909 0.005286
## ma1 0.461919 0.231952 1.99145 0.046432
## omega 0.000008 0.000032 0.26161 0.793622
## alpha1 0.069502 0.055770 1.24624 0.212676
## beta1 0.922740 0.023193 39.78587 0.000000
##
## LogLikelihood : 2857.022
##
## Information Criteria
## ------------------------------------
##
## Akaike -4.5238
## Bayes -4.4953
## Shibata -4.5239
## Hannan-Quinn -4.5131
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1083 0.7421
## Lag[2*(p+q)+(p+q)-1][8] 1.8068 1.0000
## Lag[4*(p+q)+(p+q)-1][14] 3.5229 0.9897
## d.o.f=3
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.7447 0.3882
## Lag[2*(p+q)+(p+q)-1][5] 1.7398 0.6810
## Lag[4*(p+q)+(p+q)-1][9] 2.3362 0.8613
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1378 0.500 2.000 0.7105
## ARCH Lag[5] 0.2343 1.440 1.667 0.9572
## ARCH Lag[7] 0.7626 2.315 1.543 0.9488
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 2.0886
## Individual Statistics:
## mu 0.1072
## ar1 0.1211
## ar2 0.1312
## ma1 0.1255
## omega 0.1495
## alpha1 0.1600
## beta1 0.2669
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.1037 0.9174
## Negative Sign Bias 0.4995 0.6175
## Positive Sign Bias 0.8016 0.4230
## Joint Effect 1.8672 0.6004
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 63.71 9.863e-07
## 2 30 81.10 8.000e-07
## 3 40 83.75 4.151e-05
## 4 50 98.97 3.130e-05
##
##
## Elapsed time : 0.1812401Rendimientos estandarizados
stdret <- residuals(garchfit, standardize = TRUE)chart.Histogram(stdret, methods = c("add.normal", "add.density"),
colorset = c("grey", "red", "blue"))
garchforecastt <- ugarchforecast(fitORspec = garchfit, n.ahead = 5)- Tiene sesgo positivo, significa que hay más valores positivos.
- Se tiene una distribución leptocúrtica
- Es casi simétrica
Definición del modelo utilizando la distribución t-student
garcht <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "sGARCH"),
distribution.model = "sstd")garchfitt <- ugarchfit(data = tv.L1.return, spec = garcht)
coef(garchfitt)## mu omega alpha1 beta1 skew
## -1.238879e-03 8.415623e-06 8.529305e-02 9.131079e-01 9.439485e-01
## shape
## 4.546731e+00Gráfica de la serie de rendimientos
archcoef <- coef(garchfitt)
garchunvar <- uncvariance(garchfitt)
garchmean <- fitted(garchfitt)
garchcoef## mu omega alpha1 beta1
## -1.382547e-03 8.598256e-06 7.080587e-02 9.215096e-01Ecuación: ^2_t=8.55*10^-6 + .071_t-1 $e^2 + .92^2_t-1
sqrt(garchunvar)## [1] 0.07254481garchvol <- sigma(garchfitt)
plot(garchvol)
tail(garchvol,1)## [,1]
## 2022-07-06 0.02890039Volatidades Pronosticadas
Primeras 5
sigma(garchforecastt)## 2022-07-06
## T+1 0.02733729
## T+2 0.02738946
## T+3 0.02744114
## T+4 0.02749232
## T+5 0.02754302Primeras 10
ug_forecast = ugarchforecast(garchfitt, n.ahead = 10)
ug_forecast##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 10
## Roll Steps: 0
## Out of Sample: 0
##
## 0-roll forecast [T0=2022-07-06]:
## Series Sigma
## T+1 -0.001239 0.02788
## T+2 -0.001239 0.02801
## T+3 -0.001239 0.02814
## T+4 -0.001239 0.02827
## T+5 -0.001239 0.02839
## T+6 -0.001239 0.02852
## T+7 -0.001239 0.02864
## T+8 -0.001239 0.02876
## T+9 -0.001239 0.02889
## T+10 -0.001239 0.02901Estimación de la volatilidad de la acción de Microsoft con GJR-GARCH
msftpect <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "sstd")response <- newsimpact(garchfitt)
plot(response$zx, response$zy, xlab = "prediction error", ylab = "predicted variance")
Definición modelo GARCH
garchfittnorm <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "sGARCH"),
distribution.model = "norm")Definición modelo GJR-GARCH
garchfittstd <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "sstd")Estimación de los modelos y de volatilidades. Modelo GARCH
garchfitmsft <- ugarchfit(data = tv.L1.return, spec = garchfittnorm)
garchvolmsft <- sigma(garchfitmsft)Modelo GJR-GARCH
gjrgarchfitmsft <- ugarchfit(data = tv.L1.return, spec = garchfittstd)
gjrgarchvolmsft <- sigma(gjrgarchfitmsft)Comparación de volatilidades estimadas.
plotvol <- plot(abs(tv.L1.return), col = "grey")
plotvol <- addSeries(gjrgarchvolmsft , col = "red", on=1)
plotvol <- addSeries(garchvolmsft, col = "blue", on=1)
plotvol
plotvol2 <- plot(gjrgarchvolmsft , col = "red", on=1)
plotvol2 <- addSeries(garchvolmsft, col = "blue", on=1)
plotvol2
Se realiza el pronóstico de la volatilidad para un horizonte de 10 observaciones.
garchformsft <- ugarchforecast(fitORspec = garchfitmsft, n.ahead = 10)
gjrgarchformsft <- ugarchforecast(fitORspec = gjrgarchfitmsft, n.ahead = 10)Resultados del ajuste del modelo
Coeficientes estimados.
garchcoef <- coef(garchfitmsft)
gjrgarchcoef <- coef(gjrgarchfitmsft)
garchcoef## mu omega alpha1 beta1
## -1.382547e-03 8.598256e-06 7.080587e-02 9.215096e-01gjrgarchcoef## mu omega alpha1 beta1 gamma1
## -1.347462e-03 5.444196e-06 2.884607e-02 9.339965e-01 7.415951e-02
## skew shape
## 9.460643e-01 4.762442e+00Variaza incondicional.
garchunvar <- uncvariance(garchfitmsft)
gjrgarchunvar <- uncvariance(gjrgarchfitmsft)
garchunvar## [1] 0.001118905sqrt(garchunvar)## [1] 0.03345003gjrgarchunvar## [1] 0.005443598sqrt(gjrgarchunvar)## [1] 0.07378074Volatilidades Pronosticadas
sigma(garchformsft)## 2022-07-06
## T+1 0.02733729
## T+2 0.02738946
## T+3 0.02744114
## T+4 0.02749232
## T+5 0.02754302
## T+6 0.02759323
## T+7 0.02764297
## T+8 0.02769224
## T+9 0.02774104
## T+10 0.02778938sigma(gjrgarchformsft)## 2022-07-06
## T+1 0.03018711
## T+2 0.03026209
## T+3 0.03033682
## T+4 0.03041129
## T+5 0.03048550
## T+6 0.03055946
## T+7 0.03063316
## T+8 0.03070662
## T+9 0.03077982
## T+10 0.03085278Value at Risk
La volatilidad es la medida de cuantifiación de riesgo, comunmente se realiza una estimación del VaR al 95% y al 99% en las instituciones financieras, permite estimar la perdida máxima aceptable para el 95% o 99% de los casos, este 95% o 99% coincide con el cuantíl 5% de la distribución de los retornos.
Workflow para estimación del VaR
Se especifica el tipo de modelo GARCH que utilizaremos.
gjrgarchmsft <- ugarchspec(mean.model = list(armaOrder = c(0,0)),
variance.model = list(model = "gjrGARCH"),
distribution.model = "sstd")Se especifica el tipo de ventana móvil que se utilizará
garchroll <- ugarchroll(gjrgarchmsft, data = tv.L1.return, n.start = 250,
refit.window = "moving", refit.every = 100)garchVaR <- quantile(garchroll, probs = 0.05)actual <- xts(as.data.frame(garchroll)$Realized, time(garchVaR))
VaRplot(alpha = 0.05, actual = actual, VaR = garchVaR)