Econometría Financiera
Copertino Quispe1
SUMMARY
Este es el desarrollo de la segunda evaluación del curso de Econometría Financiera realizado por EDUCATE Perú Consultores.
PREGUNTA 1
Descargar los datos del ETF de Yahoo Finance desde el 30/01/2004 hasta el 22/04/2020. El ticker del Vanguard Value ETF en Yahoo Finance es: “VTV”.
Para iniciar con la descarga del Vanguard Value ETF, se requiere cargar y llamar al paquete ; seguidamente se utiliza la función para llamar y descargar la data del , así como indica el siguiente comando:
install.packages("quantmod")
library("quantmod") getSymbols.yahoo('VTV',env=globalenv(),return.class="xts",
from='2004-01-30',to='2020-04-22',
periodicity='daily')## [1] "VTV"head(VTV)## VTV.Open VTV.High VTV.Low VTV.Close VTV.Volume VTV.Adjusted
## 2004-01-30 49.45 49.45 49.25 49.25 13200 31.95637
## 2004-02-02 49.42 49.76 49.25 49.41 2900 32.06020
## 2004-02-03 49.52 49.52 49.40 49.40 3800 32.05370
## 2004-02-04 49.30 49.30 49.03 49.19 181800 31.91744
## 2004-02-05 49.12 49.29 49.09 49.19 3100 31.91744
## 2004-02-06 49.32 49.68 49.32 49.68 7300 32.23538PREGUNTA 2
Graficar en distintas ventanas las siguientes series de tiempo:
Precio Ajustado del ETF (Adjusted Price)
Retornos del ETF
Para realizar el gráfico del Precio Ajustado del ETF, solo ploteamos el Adjusted Price del data set, para ello se requiere instalar y llamar al paquete , y utilizar la función como a continuación se verifica en el siguiente código:
install.packages("dygraph")
library("dygraph") Para realizar el gráfico de los retornos del Precio Ajustado del ETF, utilizaremos el método con la finalidad de obtener dichos retornos, para ello se requiere instalar y llamar al paquete , y utilizar la función como a continuación se verifica en el siguiente código:
install.packages("PerformanceAnalytics")
library("PerformanceAnalytics")
VTV.Ret = Return.calculate(VTV$VTV.Adjusted, method = "compound")
VTV.Ret=VTV.Ret[-1,]Finalmente realizamos el gráfico de los retornos del
PREGUNTA 3
Utilizar la estrategia de estimación para estimar la volatilidad de dicho ETF. Tener en cuenta lo siguiente:
- Describir el mejor modelo obtenido para la media condicionada de los retornos.
- Describir el mejor modelo obtenido para la varianza condicionada de los retornos.
- En caso de existir asimetría, justificar qué pruebas utilizó para detectar dicho problema y describir el mejor modelo asimétrico que usted obtuvo.
- Proyectar la varianza para 6 períodos fuera de la muestra.
Para iniciar el análisis, verificaremos si los residuos se encuentran correlacionados, para ello nos ayudaremos del gráfico de autocorrelación de los residuos, dicho gráfico se puede realizar con el siguiente código:
install.packages("PerformanceAnalytics")
library("PerformanceAnalytics") par(mfrow=c(2,2))
acf(VTV.Ret,main="ACF de los Retornos");
pacf(VTV.Ret,main="PACF de los Retornos");
acf(VTV.Ret^2,main="ACF de los Retornos^2");
pacf(VTV.Ret^2,main="PACF de los Retornos^2")
Del gráfico de las autocorrelaciones se puede mencionar que:
Existe autocorrelación debil en la ACF de los Retornos, a lo mucho un rezago, todas las demás autocorrelaciones se encuentran dentro de la banda, lo que significa que es significativamente cero, por lo que no existiría autocorrelación.
Se observa que en el gráfico de PACF de los Retornos, solo el primer rezago tiene autocorrelación fuerte negativa, y luego todas las demás se encuentran casi dentro de las bandas.
En el gráfico de ACF de los Retornos^2 se oberva que existen patrones de autocorrelación la cual cae exponencialmente.
Se oberva que en el gráfico de PACF de los Retornos^2 existen correlaciones más fuertes las cuales salen fuera de las bandas.
Luego del análisis de autocorrelación, analizaremos si los retornos se distribuyen como una normal, para ello graficaremos la distribución de los retornos y la compararemos con la función de distribución normal.
Dicho gráfico se puede realizar con el siguiente código:
m = mean(VTV.Ret) # Media
s = sd(VTV.Ret) # Desv Standar
par(mfrow=c(1,2))
hist(VTV.Ret,nclas=40,freq=F,main="Histograma de retornos");
curve(dnorm(x,mean=m,sd=s),from=-0.1,to=0.1,add=T,col="blue")
plot(density(VTV.Ret),main="Retorno de distibución empírica");
curve(dnorm(x,mean=m,sd=s),from=-0.1,to=0.1,add=T,col="blue")
Del análisis gráfico se puede apreciar que la distribución de los retornos no se comporta como una normal y que la cola del histograma de los retornos es más pesada que la Normal, ya que hay más valores extremos en la distribución del histograma de retornos que la distribución Normal.
Esto quiere decir que hay mayor probabilidad de tener valores extremos que la Normal, y que esta distribución no puede capturar dichos valores.
Lo mencionado en las líneas precedentes se puede observar en el gráfico siguiente, donde la línea roja (Normal) se encuentra por debajo del histograma de los retornos.
hist(VTV.Ret,nclas=100,freq=F,main="Histograma de retornos", xlim=c(-0.06,0),ylim=c(0,5));
curve(dnorm(x,mean=m,sd=s),from=-0.06,to=0.0,add=T,col="red")
Finalmente verificaremos el comportamiento de la disttribución de los retornos con el gráfico de QQ-Plot, donde se puede observar que las colas se están desviando mucho de la línea de 45°, ello quiere decir que la distribución no es Normal.
qqnorm(VTV.Ret,col="blue")
qqline(VTV.Ret,col="black")
A modo de conclusión, verificamos que la distribución no se comporta como una Normal por presencia de clusters de volatilidad, y dada esta característica, nos permite concluir que la volatilidad no es constante en el tiempo, por lo que se realizará la modelación de su comportamiento.
En esta sección estimaremos la volatilidad de los retornos a través del modelo EWMA (Promedio Móvil Ponderado Exponencialmente). Esta medida permite a las observaciones más recientes tener un mayor impacto (peso) en el pronóstico de la volatilidad.
El modelo EWMA puede ser expresado de varias maneras, por ejemplo:
\[ \begin{aligned} &\sigma^2_t = \lambda \sigma^2_{t-1} + (1-\lambda)(r_{t-1}-\overline{r})^2\\ &\sigma^2_t = (1-\lambda) \sum_{n=1}^{m} \lambda^{n-1}(r_{t-1}-\overline{r})^2\ \end{aligned} \]
A continuación construiremos el modelo a través de los siguientes comandos:
# i) Creamos una matriz EWMA que guarde la volatilidad (dinámica) s2
VTV.EWMA = matrix(nrow=nrow(VTV.Ret),ncol=1) # = s2
# ii) Lambda
lambda = 0.94 # Ya que estoy trabajando con datos diarios
# iii) El primer valor de la Matriz EWMA = varianza de los retornos
VTV.R2=(VTV.Ret-m)^2
VTV.EWMA[1,1]=s^2
# iv) Creamos un bucle para la ecuación EWMA que llene
# la Matriz EWMA a partir del segundo elemento
for (i in 2:nrow(VTV.EWMA)) {
VTV.EWMA[i,1]=lambda*VTV.EWMA[i-1,1]+(1-lambda)*(VTV.Ret[i-1,1]-m)^2
}
# v) Se transforma la volatilidad EWMA en formato serie de tiempo
VTV.EWMA=xts(VTV.EWMA,order.by = as.Date(index(VTV.Ret)))
names(VTV.EWMA)="EWMA"# vi) Creamos el gráfico de la volatilidad estimada EWMA (s2_EWMA)
chart.TimeSeries(VTV.EWMA,main="VTV-Modelo EWMA para la volatilidad")
Para modelar este modelo, se utilizará el estimador de rango diario como proxy de la volatilidad de los retornos.
Esta proxy tiene como forma funcional:
\[\sigma^2_t = log(\frac{alto_t}{bajo_t})\]
A continuación construiremos el modelo a través de los siguientes comandos:
# Estimador de Rango Diario
VTV.Range=log(VTV$VTV.High/VTV$VTV.Low)
names(VTV.Range)="Range"chart.TimeSeries(VTV.Range,main="VTV-Estimador de rango diario")
Finalmente concluímos que el modelo EWMA es más suavizado que el de Rango de Volatilidad; a su vez, el Modelo EWMA es quien captura mejor las volatilidades recientes producto de la pandemia del coronavirus, efecto que no es capturado por el estimador de rango diario
Para modelar estos modelos, se utilizará la estrategia de estimación con la finalidad de encontrar el mejor modelo para la varianza condicionada de los retornos.
Paso 1)
Encontrar el mejor modelo ARMA(p,q) (via Box-Jenkins), verificando que el error es ruido blanco.
Como ya se había verificado en el gráfico de las autocorrelaciones de los retornos, encontramos que en su ACF a lo mucho se cuenta con un rezago y en su PACF solo el primer rezago tiene autocorrelación fuerte .
También debemos recordar que según la gráfica de los retornos, esta es estacionaria en media (= 0), por lo que no se diferenciaria ni se le pasará el test de Raíz Unitaria.
Estos resultados servirán como insumo para encontrar el mejor modelo ARMA(p,q).
Para iniciar con la modelación, se requiere cargar y llamar al paquete ; seguidamente se utiliza la función y se completa los argumentos ya analizados, así como indica el siguiente comando:
install.packages("tseries")
library("tseries") A continuación construiremos el modelo con los siguientes comandos:
VTV.ARIMA=auto.arima(VTV.Ret,d=NA,D=NA,max.d=0,max.D=NA,
max.p=1,max.q=1,max.P=NA,max.Q=NA,
seasonal=F,xreg=NULL,ic=c("aic"))
VTV.ARIMA## Series: VTV.Ret
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 mean
## -0.1253 3e-04
## s.e. 0.0153 2e-04
##
## sigma^2 estimated as 0.0001538: log likelihood=12131.12
## AIC=-24256.25 AICc=-24256.24 BIC=-24237.3Este resultado nos muestra que el mejor modelo estimado para la media condicionada es el .
Paso 2)
Verificar Clusters de volatilidad.
Para verificar el comportamiento de los residuos los plotaremos y analizaremos el test de Ljung Box
Para iniciar con la modelación, se requiere cargar y llamar al paquete ; seguidamente se utiliza la función , así como indica el siguiente comando:
install.packages("forecast")
library("forecast") checkresiduals(VTV.ARIMA) 
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,1) with non-zero mean
## Q* = 42.992, df = 8, p-value = 8.815e-07
##
## Model df: 2. Total lags used: 10Según el Test de Ljung Box se acepta la Ho de que no hay correlación serial.
Asimismo, según el correlograma de los residuos al cuadrado, se puede observar que las barras de la autocorrelación simple no sobresalen de las bandas (exceptuando el lag 4); asi también, en la autocorrelación parcial se puede observar que no hay correlación serial (o bueno si la hay es en el orden 4 o de mayor orden), pero en las ordenes iniciales no se visualiza ello. Concluyendo que no hay existe correlación serial de los residuos.
par (mfrow=c(1,2))
acf (VTV.ARIMA$residuals,main="ACF de los Residuos");
pacf(VTV.ARIMA$residuals,main="PACF de los Residuos")
Ahora verificaremos el comportamiento de los cluster de volatilidad, ploteando los residuos al cuadrado.
par (mfrow=c(1,2))
acf (VTV.ARIMA$residuals^2,main="ACF de los Residuos^2");
pacf(VTV.ARIMA$residuals^2,main="PACF de los Residuos^2")
Concluimos que sí hay Cluster de volatilidad, ya que las funciones de autocorrelaciones no se eliminan (ni la acf ni pacf), verificando la existencia de patrones para la varianza, la cual se modelará a través del modelo .
Paso 3)
Una vez identificado los Clusters de volatilidad, seleccionaremos los modelos GARCH candidatos, luego de observar el correlograma de los residuos al cuadrado.
Para ello, se programará una rutina con la finalidad de obtener el mejor modelo GARCH candidato.
# Programa para elegir el mejor modelo GARCH
GARCH.Mods = list()
GARCH.Fit = list()
GARCH.BIC = list()
# Programamos un bucle para realizar los modelos
# Recordemos del autoarima p.max=1 y q.max=1
for (q in 1:2) {
for (p in 1:2) {
GARCH.Mods[[paste("GARCH",q,p,sep="")]]=
ugarchspec(variance.model=list(model="sGARCH",
garchOrder=c(q,p)),mean.model = list(armaOrder=c(0,1)))
GARCH.Fit[[paste("GARCH",q,p,sep="")]]=
ugarchfit(spec=GARCH.Mods[[paste("GARCH",q,p,sep="")]],data=VTV.Ret)
GARCH.BIC[[paste("GARCH",q,p,sep="")]]=
infocriteria(GARCH.Fit[[paste("GARCH",q,p,sep="")]])
GARCH.BIC[[paste("GARCH",q,p,sep="")]]=
GARCH.BIC[[paste("GARCH",q,p,sep="")]][2,1]
}
}
# Data Frame de los modelos
GARCH.BIC=as.data.frame(GARCH.BIC)
min.GARCH=min(GARCH.BIC)
bic.GARCH=match(min.GARCH,GARCH.BIC)
# El mejor modelo GARCH seleccionado:
colnames(GARCH.BIC[match(min(GARCH.BIC),GARCH.BIC)])## [1] "GARCH11"Finalmente, del análisis previo, concluimos que el mejor modelo es el .
Paso 4)
Una vez evaluado el mejor modelo GARCH, verificaremos que los errores estandarizados estén limpios de volatilidad.
Para ello evaluaremos el correlograma del los errores estandarizados a través de los siguientes comandos:
par(mfrow=c(2,2))
acf (residuals(GARCH.Fit$GARCH11,standardize=T),main="ACF de los Residuos")
pacf(residuals(GARCH.Fit$GARCH11,standardize=T),main="ACF de los Residuos")
acf (residuals(GARCH.Fit$GARCH11,standardize=T)^2,main="ACF de los Residuos^2")
pacf(residuals(GARCH.Fit$GARCH11,standardize=T)^2,main="ACF de los Residuos^2")
Análisis: Practicamente ya no existen patrones de volatilidad (“está limpio”) ya que la mayoría de las correlaciones están dentro de las bandas, finalmente realizaremos realizaremos el Test de Ljung Box para contrastar las correlaciones, del cual concluimos que los residuos están limpios de cluster de volatilidad.
Box.test(residuals(GARCH.Fit$GARCH11,standardize=T),lag=10,type=c("Ljung-Box"),fitdf=0)##
## Box-Ljung test
##
## data: residuals(GARCH.Fit$GARCH11, standardize = T)
## X-squared = 17.558, df = 10, p-value = 0.06289Paso 5)
Explorar la existencia de asimetría en la NIC (Curva de impacto de las noticias).
Al evaluar la , esta no es significativa, por lo que no debería haber asimetría.
GARCH.Fit$GARCH11##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000707 0.000108 6.5213 0.000000
## ma1 -0.039291 0.017475 -2.2484 0.024550
## omega 0.000002 0.000001 3.0503 0.002286
## alpha1 0.140397 0.011470 12.2398 0.000000
## beta1 0.840068 0.011908 70.5463 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000707 0.000101 7.01455 0.000000
## ma1 -0.039291 0.016183 -2.42796 0.015184
## omega 0.000002 0.000005 0.52316 0.600862
## alpha1 0.140397 0.026590 5.27997 0.000000
## beta1 0.840068 0.045408 18.50056 0.000000
##
## LogLikelihood : 13554.9
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.6372
## Bayes -6.6295
## Shibata -6.6372
## Hannan-Quinn -6.6345
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5257 0.4684
## Lag[2*(p+q)+(p+q)-1][2] 0.5332 0.9559
## Lag[4*(p+q)+(p+q)-1][5] 2.3845 0.5995
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.006531 0.9356
## Lag[2*(p+q)+(p+q)-1][5] 1.872381 0.6489
## Lag[4*(p+q)+(p+q)-1][9] 2.517573 0.8351
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.000435 0.500 2.000 0.9834
## ARCH Lag[5] 0.515361 1.440 1.667 0.8791
## ARCH Lag[7] 0.762802 2.315 1.543 0.9488
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 16.6122
## Individual Statistics:
## mu 0.03327
## ma1 0.08906
## omega 1.38626
## alpha1 0.34477
## beta1 0.67921
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.3574 0.174734
## Negative Sign Bias 0.2572 0.797027
## Positive Sign Bias 2.2984 0.021588 **
## Joint Effect 15.4194 0.001491 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 153.4 4.768e-23
## 2 30 176.2 5.151e-23
## 3 40 194.6 1.499e-22
## 4 50 213.7 1.874e-22
##
##
## Elapsed time : 0.4810278Asimismo, graficaremos la curva de impacto de las noticias para evaluar su asimetría. Para realizar dicho gráfico se requiere cargar y llamar al paquete ; seguidamente se utiliza la función para obetener las NIC del modelo estimado , así como indica el siguiente comando:
install.packages("rugarch")
library("newsimpact") NIC.GARCH11=newsimpact(GARCH.Fit$GARCH11)
plot(NIC.GARCH11$zx,NIC.GARCH11$zy,
type="l",lwd=2,col="blue",
main="GARCH(1,1)-News Impact",
ylab=NIC.GARCH11$yexpr,xlab=NIC.GARCH11$Xexpr)
Graficamente se puede observar que la Curva de Impacto de Noticias es simétrica por lo que debemos cambiar el modelo por un EGARCH o TGARCH.
Paso 6)
Si la NIC es simétrica cambiar el GARCH por un EGARCH o un TGARCH
Ya que evaluamos el Test de Signos y el grafico de la NIC concluimos que es simétrica, por lo que debemos cambiar el GARCH a un o un .
Para realizar dichos modelos, se utilizará el siguiente comando:
# Modelo EGARCH
EGARCH.Spec=ugarchspec(variance.model=list(model="eGARCH",garchOrder=c(1,1))
,mean.model = list(armaOrder=c(0,1)))
EGARCH.Fit=ugarchfit(EGARCH.Spec,data=VTV.Ret)
# Modelo TGARCH
TGARCH.Spec=ugarchspec(variance.model=list(model="fGARCH",submodel="TGARCH",
garchOrder=c(1,1)),mean.model = list(armaOrder=c(0,1)))
TGARCH.Fit=ugarchfit(TGARCH.Spec,data=VTV.Ret)Finalmente, graficaremos la curva de impacto de las noticias de los modelos GARCH. Para realizar dicho gráfico se seguirá el siguiente comando:
NIC.EGARCH11=newsimpact(EGARCH.Fit)
NIC.TGARCH11=newsimpact(TGARCH.Fit)# GARCH11
par(mfrow=c(1,1))
plot(NIC.GARCH11$zx,NIC.GARCH11$zy,type="l",lwd="2",col="blue",
ylab=NIC.GARCH11$yexpr,xlab=NIC.GARCH11$xexpr,ylim=c(0,0.009))
par(new=T)
# EGARCH11
plot(NIC.EGARCH11$zx,NIC.EGARCH11$zy,type="l",lwd="2",col="red",
ylab="",xlab="",ylim=c(0,0.009))
par(new=T)
# TGARCH11
plot(NIC.TGARCH11$zx,NIC.TGARCH11$zy,type="l",lwd="2",col="green2",
ylab="",xlab="",ylim=c(0,0.009))
title(main="News Impact Curce",font.main=4)
grid()
box()
legend(x=0.17,y=0.0035,legend=c("GARCH","EGARCH","TGARCH"),
fill=c("blue","red","green2"),cex=0.9)
Análisis: Los modelos EGARCH y TGARCH capturan mejor la variabilidad de los retornos, en el gráfico se observa que estas NIC capturan mejor los efectos asimetricos que el modelo GARCH.
TGARCH
Paso 7)
Ahora que ya capturamos los efectos asimetricos de los choques verificaremos la Normalidad de los residuos de estos nuevos modelos asimetricos.
Para ello, cargaremos y llamaremos al paquete ; seguidamente se utiliza la función para analizar el test de Jarque Bera sobre la normalidad de los residuos, así como indica el siguiente comando:
install.packages("normtest")
library("jb.norm.test") jb.norm.test(residuals(TGARCH.Fit,standardize=T))##
## Jarque-Bera test for normality
##
## data: residuals(TGARCH.Fit, standardize = T)
## JB = 912.29, p-value < 2.2e-16Según el test de JB (p-value < 2.2e-16) rechazamos la Ho de que los residuos estandarizados no se comportan como una Normal.
Paso 8)
Si el error estandarizado no es normal, entonces cambiar a distribucion t o GED.
Dado que los residuos no son normales, cambiaremos la especificación del modelo por una distribución T-Student (std).
TGARCH.Spec=ugarchspec(variance.model=
list(model="fGARCH",submodel="TGARCH",garchOrder=c(1,1)),
mean.model = list(armaOrder=c(0,1)),distribution.model="std") Cuando especificamos distribution.model=“std”, vemos al inicio del modelo estimado que la distribución ahora es t-studen (std)
TGARCH.Fit=ugarchfit(TGARCH.Spec,data=VTV.Ret)
TGARCH.Fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : fGARCH(1,1)
## fGARCH Sub-Model : TGARCH
## Mean Model : ARFIMA(0,0,1)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000524 0.000102 5.1184 0.000000
## ma1 -0.036891 0.015470 -2.3847 0.017093
## omega 0.000232 0.000033 6.9769 0.000000
## alpha1 0.101601 0.010200 9.9604 0.000000
## beta1 0.897402 0.009438 95.0852 0.000000
## eta11 0.979749 0.100178 9.7801 0.000000
## shape 6.426509 0.632484 10.1607 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000524 0.000100 5.2400 0.00000
## ma1 -0.036891 0.013725 -2.6879 0.00719
## omega 0.000232 0.000044 5.3353 0.00000
## alpha1 0.101601 0.012354 8.2242 0.00000
## beta1 0.897402 0.012035 74.5671 0.00000
## eta11 0.979749 0.109005 8.9881 0.00000
## shape 6.426509 0.645569 9.9548 0.00000
##
## LogLikelihood : 13748.96
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.7313
## Bayes -6.7205
## Shibata -6.7313
## Hannan-Quinn -6.7275
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1600 0.6891
## Lag[2*(p+q)+(p+q)-1][2] 0.5009 0.9641
## Lag[4*(p+q)+(p+q)-1][5] 2.0315 0.7022
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.08914 0.7653
## Lag[2*(p+q)+(p+q)-1][5] 1.51080 0.7372
## Lag[4*(p+q)+(p+q)-1][9] 2.17382 0.8833
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.001192 0.500 2.000 0.9725
## ARCH Lag[5] 0.471222 1.440 1.667 0.8921
## ARCH Lag[7] 0.814357 2.315 1.543 0.9418
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 3.694
## Individual Statistics:
## mu 0.39165
## ma1 0.04544
## omega 0.76229
## alpha1 0.89344
## beta1 1.36749
## eta11 0.73098
## shape 0.68667
##
## 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.5521 0.58091
## Negative Sign Bias 1.0134 0.31094
## Positive Sign Bias 1.9420 0.05221 *
## Joint Effect 4.8210 0.18539
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 75.71 1.008e-08
## 2 30 79.92 1.191e-06
## 3 40 100.21 2.739e-07
## 4 50 112.21 7.277e-07
##
##
## Elapsed time : 1.458084EGARCH
Paso 7)
Ahora verificaremos la Normalidad de los residuos para el modelo EGARCH.
jb.norm.test(residuals(EGARCH.Fit,standardize=T))##
## Jarque-Bera test for normality
##
## data: residuals(EGARCH.Fit, standardize = T)
## JB = 996.11, p-value < 2.2e-16Según el test de JB (p-value < 2.2e-16) rechazamos la Ho de que los residuos estandarizados, concluyendo que no se comportan como una Normal.
Paso 8)
Si el error estandarizado no es normal, entonces cambiar a distribucion t o GED
Dado que los residuos no son normales, cambiaremos la especificación del modelo por una distribución T-Student (std).
EGARCH.Spec=ugarchspec(variance.model=
list(model="eGARCH",garchOrder=c(1,1)),
mean.model = list(armaOrder=c(0,1)),distribution.model="std") Cuando especificamos distribution.model=“std”, vemos al inicio del modelo estimado que la distribución ahora es t-studen (std)
EGARCH.Fit=ugarchfit(EGARCH.Spec,data=VTV.Ret)
EGARCH.Fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : eGARCH(1,1)
## Mean Model : ARFIMA(0,0,1)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000592 0.000093 6.3656 0.000000
## ma1 -0.039345 0.015692 -2.5073 0.012166
## omega -0.206594 0.002944 -70.1637 0.000000
## alpha1 -0.159590 0.008798 -18.1404 0.000000
## beta1 0.978905 0.000094 10406.8562 0.000000
## gamma1 0.176314 0.010290 17.1347 0.000000
## shape 6.158806 0.586166 10.5069 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.000592 0.000091 6.4933 0.000000
## ma1 -0.039345 0.013609 -2.8912 0.003838
## omega -0.206594 0.003340 -61.8505 0.000000
## alpha1 -0.159590 0.010238 -15.5875 0.000000
## beta1 0.978905 0.000109 8943.1000 0.000000
## gamma1 0.176314 0.011785 14.9611 0.000000
## shape 6.158806 0.621700 9.9064 0.000000
##
## LogLikelihood : 13734.9
##
## Information Criteria
## ------------------------------------
##
## Akaike -6.7244
## Bayes -6.7136
## Shibata -6.7244
## Hannan-Quinn -6.7206
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.2271 0.6337
## Lag[2*(p+q)+(p+q)-1][2] 0.5654 0.9467
## Lag[4*(p+q)+(p+q)-1][5] 2.2748 0.6314
## d.o.f=1
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.02595 0.8720
## Lag[2*(p+q)+(p+q)-1][5] 1.29087 0.7911
## Lag[4*(p+q)+(p+q)-1][9] 1.79682 0.9278
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.03849 0.500 2.000 0.8445
## ARCH Lag[5] 0.38843 1.440 1.667 0.9159
## ARCH Lag[7] 0.54631 2.315 1.543 0.9740
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 3.6168
## Individual Statistics:
## mu 0.20628
## ma1 0.04928
## omega 0.88824
## alpha1 0.20388
## beta1 0.71969
## gamma1 0.77036
## shape 0.37271
##
## 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.2732 0.78473
## Negative Sign Bias 0.6634 0.50711
## Positive Sign Bias 1.7658 0.07751 *
## Joint Effect 3.6416 0.30286
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 62.67 1.451e-06
## 2 30 68.18 5.365e-05
## 3 40 89.71 7.187e-06
## 4 50 108.71 2.033e-06
##
##
## Elapsed time : 1.025059Finalmente agruparemos la volatilidad estimada por cada modelo estimado.
GARCH.Vol=cbind(GARCH.Fit$GARCH11@fit$sigma,TGARCH.Fit@fit$sigma,EGARCH.Fit@fit$sigma)
GARCH.Vol=xts(GARCH.Vol,order.by = index(VTV.Ret))
names(GARCH.Vol)=c("GARCH","TGARCH","EGARCH")Y graficaremos lo resultados, para ello cargaremos y llamaremos al paquete ; seguidamente utilizaremos la función para realizar el gráfico combinado, así como indica el siguiente comando:
install.packages("PerformanceAnalytics")
library("chart.TimeSeries") chart.TimeSeries(GARCH.Vol,main="Volatilidades estimadas",legend.loc="topleft",lwd=1)
Concluimos que la curva de impactos de noticias de los modelos y capturan mejor los efectos asimetricos de un shock sobre la volatilidad.
Además, se observa que ante un shock negativo existe mayor incertidumbre por lo que la volatilidad es mayor, a diferencia de un shock positivo, esto para los modelos y , ya que para el modelo los choques negativos y positivos tienen el mismo efecto sobre la volatilidad.
Por ello los modelos TGARCH y EGARCH, capturan mejor las volatilidades del retorno.
A modo de conclusión, es el modelo TGARCH el que mejor captura las volatilidades del modelo como se puede observar en el gráfico, por lo que nos quedaremos con el modelo como el modelo estimado.
Paso 9)
Concluiremos el siguiente informe con la proyección de 6 períodos hacia adelante con el modelo estimado .
A continuación presentamos los códigos para su realización:
# Proyección de 6 períodos
GARCH.Forecast = ugarchforecast(GARCH.Fit$GARCH11, data = SPY.Ret, n.ahead = 6)
TGARCH.Forecast = ugarchforecast(TGARCH.Fit, data = SPY.Ret, n.ahead = 6)
EGARCH.Forecast = ugarchforecast(EGARCH.Fit, data = SPY.Ret, n.ahead = 6)
VOL.Forecast = cbind(GARCH.Forecast@forecast$sigmaFor,
TGARCH.Forecast@forecast$sigmaFor,
EGARCH.Forecast@forecast$sigmaFor)
VOL.Forecast = ts(VOL.Forecast)
# Gráfico de proyecciones
par(mfrow = c(1,1))
plot(VOL.Forecast[,1], type="l", lwd=2, col = "blue",
ylab=bquote(sigma[t]^2), ylim = c(0.02,0.055))
par(new = TRUE)
plot(VOL.Forecast[,2], type = "l", lwd = 2, col = "red",
ylab=bquote(sigma[t]^2), ylim = c(0.02,0.055))
par(new = TRUE)
plot(VOL.Forecast[,3], type = "l", lwd = 2, col = "green2",
ylab=bquote(sigma[t]^2), ylim = c(0.02,0.055))
title(main = bquote(sigma[t]^2 ~ "forecast"), font.main = 4)
grid()
box()
legend(x = 5, y = 0.05, legend = c("GARCH", "TGARCH", "EGARCH"),
fill = c("blue", "red", "green2"), cex = 0.9)
UNIVERSIDAD NACIONAL MAYOR DE SAN MARCOS, coperquispe@gmail.com↩