Bienvenidos al tutorial de Modelos Heterocedásticos y Series de Rendimientos. Los objetivos de este tutorial es que aplicaremos conceptos de los modelos heterocedásticos, particularmente derivaremos un modelo GARCH (p,q) (1,1). Para ello, vamos a obtener la serie de precios de cierre del activo SPY, que es un fondo indexado al estándar Plus500 entre las fechas 2017 y 2021, y vamos a obtener sus rendimientos simples. Vamos a generar un modelo ARMA(p,q) inicialmente para ver cuál es ese modelo de media de rendimientos y vamos a revisar los residuales para probar si hay, efectivamente, el efecto ARCH. Como si lo va haber, vamos a generar este modelo GARCH(1,1) con una especificación particular a ese comportamiento de rendimientos, en cuanto a su sesgo y su curtosis. Vamos a revisar si es un modelo con reversión a la media y con qué valor, vamos a generar un pronóstico y, por último, una simulación de sus parámetros para las siguientes fechas. Lo que vamos a obtener, particularmente, es que vamos a obtener una serie rendimientos con sesgo significativo y un exceso de curtosis positiva significativa. Esta se va a ver de esta manera así. Como ustedes pueden ver, es una distribución leptocúrtica con presencia de colas anchas y, en este caso, también vamos a tener un sesgo significativo. Como esta distribución va a responder a una distribución de colas pesadas y sesgadas, la vamos a poder modelar, como tal, con una distribución en las innovaciones como una t-student sesgada. Y es lo que vamos a alimentar en el modelo de volatilidad. Como podemos ver, esta es una distribución normal, la línea roja, la distribución de referencia con su misma media y varianza. Esta sería una distribución normal con esos parámetros; esta es la distribución de rendimientos y esto, en gris, su histograma. ¿Cómo funciona la modulación? Vamos a emplear el modelo rugarch en R, donde requiere de tres pasos. Los primeros dos a continuación. El primero es donde necesitamos dar una especificación al modelo, en tanto la media, la varianza y la distribución de las innovaciones o de los choques, que comentábamos se veían en esa distribución leptocúrtica. El segundo paso es estimar el modelo, el ajuste Fit, con esa especificación inicial. De esta manera, si nosotros vemos el código, lo primero que tenemos que hacer es alimentar el modelo de media, que vamos a encontrar en una primera parte. Para ello, vamos a ver que es un modelo ma2. Lo podemos modelar aquí y decir que sí incluimos la media, si esta media es diferente o igual a cero. Lo segundo, especificar el modelo de varianza. En este caso, vamos a trabajar con un estándar GARCH y vamos a trabajar con un GARCH(1,1). Lo segundo es que realmente el primer paso de este rugarch es alimentar estos dos pasos, estas dos definiciones de modelo de media y el modelo de la varianza, en esta UGARCH "especification", "ugarchspec", como dice aquí. Si usted pregunta, "¿Cuál es ese modelo de varianza?", vamos a agarrar el modelVar1, que ya dijimos arriba cuál iba a ser. Igualmente, el modelo de media que definimos, que era un ma2 y, en este caso particular, dado este comportamiento sesgado y de distribuciones con colas anchas, podemos elegir una distribución t sesgada. Esa se llama "sstd". Si no hubiese sido sesgada, únicamente con distribuciones de colas anchas, hubiera sido "std". El segundo paso es el modelFit y podemos generar el ajuste. Si nosotros solicitamos los coeficientes a ello, vamos a tener mu, ma1, ma2, que es el modelo de media, omega, alfa y beta, y estos dos parámetros que se refieren a la distribución t sesgada. Como tenemos un modelo ARMA(0,2) y un modelo GARCH(1,1) con innovaciones t student sesgados, lo podemos modelar a partir de estos resultados de los parámetros de esta manera. El módulo de media o la ecuación de media de rendimiento está dado por estos parámetros que nos acaba de entregar, y el modelo de volatilidad está dado por un GARCH(1,1), donde tenemos los coeficientes alfa, beta y los dos parámetros para la distribución t sesgada. El tercer paso sería la previsión a partir de este modelo de ajuste. Podemos solicitar la predicción a tantos pasos adelante y esta sería una parte de la serie, y cuál sería ese pronóstico con estas bandas de intervalos de confianza. Este es el actual y este el forecast. Con respecto particular al comportamiento de la sigma, tendríamos la desviación estándar, el comportamiento actual y el pronóstico de un pedazo de esas series que tenemos por aquí. Si quisiéramos hacer esa simulación, podemos agarrar a partir de dos parámetros y hacer una simulación de cuál sería el comportamiento de nuestra serie, tanto para los valores de la varianza condicional y su distribución. Cómo se verían estos dos parámetros, este parámetro de varianza condicional en cuanto a su puro valor en el tiempo y su distribución. Como podemos ver aquí, es que esta tiene una reversión a la media, con un cierto valor sigma que vamos a determinar el código a largo plazo. Está teniendo esta reversión a la media con este tipo de distribuciones y, en cuanto a la serie de rendimientos, este es el valor actual y cómo se vería a futuro esta serie de rendimientos. Igualmente, la distribución de esos rendimientos. Pasemos al código para encontrar estos resultados. Ya nos encontramos en RStudio. Lo primero que tenemos que hacer es asegurarnos que tenemos las paqueterías instaladas y activas. Como siempre, estoy llamando a que no deshabiliten los "warnings" para no llenar la pantalla. Tenemos la librería rugarch, que es la más importante en este caso, la de pronóstico, la de TSstudio para hacer la prueba de "extended autocorrelation function" y derivar el módulo de media, el quantmod para obtener los datos de internet y demás. Esto es para encontrar el efecto ARCH, como vamos a ver más adelante. El primer paso sería la obtención de los datos y el modelo de medias. Si ustedes se dan cuenta, es con las fechas solicitadas. Esta función la estamos llamando "rendimientos", que es la misma que trabajamos en el módulo de "momentos". Esta no tiene mayor complicación. Estamos trabajando nuevamente con los datos de precio de cierre y estamos calculando los rendimientos simples. Aunque es un solo activo, si ustedes tuvieran muchos más, aquí los pueden llamar y esta función se facilita. Ya que la tenemos hecha, la vamos a llamar para nuestro simple activo del fondo indexado, el estándar Plus500, y podemos ver que se trata de un formato igualmente de xts, y aquí nos vienen los demás detalles de desde dónde, hasta qué fechas y cuántos valores está considerando. Podemos generar también un gráfico que nos dé idea si existen "clusters" de volatilidad. Como podemos ver, aquí tenemos un "cluster", aquí tenemos otro "cluster" de volatilidad a lo largo de la serie de rendimientos. Podemos pasar, con esto en mente, a modelar esos rendimientos de la serie. Vamos a aplicar metodología de Extended Autocorrection Function y auto arima para ver si tenemos los mismos resultados. Para el "extended", necesitamos la librería TSA y Basics. Para ello, necesitamos pasar nuestros datos a un DataFrame, y especificar que queremos la columna de esos rendimientos. Y podemos trabajar y obtener este. ¿Qué es lo que buscamos? Sería una esquina de un triángulo o un punto. En este caso, observamos que aquí arriba tenemos un buen candidato, que sería un ma2. Si esto fuese un tache, una cruz, a lo mejor podríamos agarrar este otro, pero este aparenta ser mejor. Esto sería un ma2 con el Extended Autocorrection Function y, si nosotros trabajamos con un auto arima, también nos estaría entregando el mismo resultado con un "non zero mean", con un resultado diferente, con una media diferente de cero, y nos está entregando sus métricos. De esta manera, tenemos los mismos resultados. Si hubieran ustedes tenido diferentes, aquí los pudieran haber modelado. Lo que podemos observar es que no tenemos problemas en los residuales. Sin embargo, podemos tener un efecto heterocedástico y esto no lo vamos a observar aquí. Vamos a llamar a esta función, a la librería aTSA y vamos a correr la función arch.test para ver si tenemos una dependencia entre los rezagos de los residuos, obviamente. Vamos a jalar este, el arch.testm y vemos que es una prueba de heterocedasticidad a los residuales, donde la hipótesis alternativa es que son heterocedásticos. Quiere decir que se puede correr a través de la prueba Portmanteau o la prueba de Langrage Multiplier, y lo que estamos obteniendo es que en ambos casos para los rezagos 4, 8, 12, 16, 20, rezagos bajos y altos, estamos teniendo un efecto ARCH, un efecto heterocedástico. Lo que podemos ver, a través de estos resultados estadísticos de prueba de Portmanteau o de Lagrange, es que tenemos "p values" muy pequeños, lo que estaría rechazando la hipótesis nula y estaría saltando a la alternativa que tenemos efectos heterocedásticos en los rezagos 4, 8, etcétera, hasta el 24. Es decir, en rezagos de residuales bajos y altos. Con eso en mente, si hubiésemos tenido rezagos muy bajos, como un 4 y los demás estos sin problema, pudiéramos haber pensado en modelarlo con un efecto ARCH. Sin embargo, como tenemos rezagos altos, significativos, una buena resolución es agarrar modelos GARCH, donde entre sus componentes estamos integrando pronósticos anteriores de la volatilidad. Con eso en mente, podemos pasar a las estadísticas básicas para darnos idea de qué modelo emplear para las innovaciones. Con eso en mente, podemos solicitar la distribución original, que tenemos una distribución sesgada y hubo un exceso de curtosis positiva. Vamos a identificar en el caso del sesgo, que pudiera haber duda por la magnitud, si es estadísticamente significativo o no, y podemos analizar la prueba t para ver si es significativamente diferente de cero. Y así es. Eso quiere decir que el sesgo es significativo. Y tenemos un exceso de curtosis, por eso es que estamos proponiendo para modelar esa serie una distribución particular, como dijimos en la presentación. Si generamos ese gráfico de densidad a sus distribuciones, obtenemos lo que estábamos viendo en la presentación. Obtenemos una serie de rendimientos leptocúrtica con sesgo negativo y exceso de curtosis positiva. Como decimos aquí, contamos con distribución de colas pesadas, sesgo significativo, podemos emplear una distribución t sesgada que requerirá dos parámetros adicionales a la distribución normal, que es shape y épsilon. Para este, el shape, a valor más bajo, quiere decir que más gordas van a ser las colas, y épsilon significa o da más detalle a la magnitud del sesgo de la distribución. En el caso de épsilon, la distribución es t sesgada, que se va a denominar en este paquete "Sstd". Si la épsilon es igual uno, es que es perfectamente simétrico. Si el épsilon es mayor a uno, es que tenemos sesgo positivo. Y, si épsilon es menor a uno, es que tenemos sesgo negativo. Nosotros tenemos un sesgo negativo, y es lo que vamos a obtener seguramente en los parámetros cuando modelemos ese ajuste a los datos. Vamos a hacer ese modelamiento como una tercera parte de este código en el modelo de volatilidad. Recordar que para un modelo de volatilidad empleando rugarch, necesitamos tres cosas. Primero es hacer la especificación, donde decimos qué modelo para la media y para la varianza tenemos. Es el ajuste dado esa especificación y el forecast. Esto lo habíamos comentado en la presentación. Corremos este segmento, que ya lo habíamos mostrado. Vamos a correrlo y nos entrega los parámetros. Lo que podemos observar a través del resultado de ajuste de esta es lo siguiente. Nos entrega cuáles son esos parámetros, tanto de la media como del modelo de volatilidad. Nos está diciendo, por ejemplo, que este ma2, este parámetro no es significativo, al igual que el valor del intercepto de ese modelo de volatilidad. ¿Qué más? Seguimos avanzando. En cuanto a los resultados ya de ese modelo y cómo están modelando los residuales, estamos teniendo que tenemos resultados de "Standardized Residuals" y "Standardized Residuals cuadráticos". ¿Qué es lo importante aquí? Que la hipótesis nula nos está diciendo que no hay una correlación serial ya de esos residuos. ¿Qué estamos teniendo? Estamos teniendo un problema, porque nos está diciendo que ese p-valor al primer rezago, al quinto y al nueve están siendo significativos. Estamos rechazando esta hipótesis nula con este valor de p-valor y estamos diciendo que sí hay correlación serial, lo cual es un problema. ¿Qué es una solución? Regresarnos a este modelo y nosotros propusimos que era una distribución t sesgada, podemos modelarlo como una distribución t no sesgada y vamos a ver qué pasa. En este sentido, vamos a correrlo. Y vamos a ver si ese modelo mejora. Si nosotros lo corremos con un modelo de distribución en las innovaciones de std, estamos obteniendo un mejor modelamiento, porque lo que está entregando es que ya no hay una correlación serial en los residuos. Con estos valores de p-valor mayores a alfa, no estamos rechazando esta hipótesis. En tanto a estas pruebas de Ljung Box en los residuos estandarizados, en los residuos estandarizados cuadráticos. Por lo tanto, una manera mejor de modular esto sería con este tipo de distribución std. Si nosotros continuamos, seguimos avanzando, podemos tener una representación de esa media de rendimientos, ahora con este modelo, y esa volatilidad, lo cual está mostrándonos que tenemos una regresión a la media de esa volatilidad. Esa volatilidad a largo plazo la podemos también calcular y este nos va a dar un 5 por ciento, pero el valor de la volatilidad es 0,002, que ese sería este valor, en cuanto está teniendo a su media, y su desviación estándar a largo plazo sería del 5 por ciento de los rendimientos. En cuanto a la simulación, tendríamos esa representación de esas simulaciones a largo plazo, como comentamos en la presentación. A partir de los resultados que obtuvimos del modelo de volatilidad con innovaciones t-student no sesgadas, obtendríamos esto, donde la diferencia notoria es aquí, en cuanto a esa distribución del error con los parámetros de la distribución t. De la misma manera, tomaríamos los mismos parámetros que obtuvimos, los modelaríamos. Esto sería un modelo ARMA(0,2), un ma2, con GARCH(1,1), con distribución de innovaciones t-student.
