UNIVERSIDAD DE COSTA RICA SISTEMA DE ESTUDIOS DE POSGRADO ESTIMACIÓN DE PARÁMETROS EN UN MODELO DE VOLATILIDAD ESTOCÁSTICA CON MEMORIA LARGA USANDO FILTRO DE PARTÍCULAS Tesis sometida a la consideración de la Comisión del Programa de Estudios de Posgrado en Matemática para optar al grado y t́ıtulo de Maestŕıa Académica en Matemática Aplicada. ANDRÉS QUIRÓS GRANADOS Ciudad Universitaria Rodrigo Facio, Costa Rica 2021 DEDICATORIA Dedicado a mi hijo Isaac, mi esposa Jennifer y mi madre Gladys. i AGRADECIMIENTOS Agradezco a mi Dios por permitirme este gran logro. “Porque Jehová da la sabi- duŕıa, y de su boca viene el conocimiento y la inteligencia”. Proverbios 2:6. Un sincero agradecimiento a mi tutor Dr. Luis Barboza por su orientación y por sus consejos los cuales ayudaron a cumplir con los objetivos. Agradezco a mi esposa por todo el apoyo incondicional, por su motivación para poder terminar el proyecto y por su ayuda siendo una lectora más del trabajo. ii TABLA DE CONTENIDOS Dedicatoria i Agradecimientos ii Hoja de aprobación iii Resumen vii Lista de cuadros viii Lista de figuras xi Lista de algoritmos xiv 1. Introducción 1 1.1. Antecedentes y justificación . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1. Objetivo general . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2. Objetivos espećıficos . . . . . . . . . . . . . . . . . . . . . . . . 4 iv 2. Marco teórico 5 2.1. Proceso estocástico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1. Proceso estocástico de memoria larga . . . . . . . . . . . . . . . 7 2.1.2. Movimiento browniano fraccionario y ruido gaussiano fraccionario 8 2.2. Modelo de volatilidad estocástica . . . . . . . . . . . . . . . . . . . . . 9 2.3. Modelo espacio estado . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4. Filtro de part́ıculas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4.1. Estimación de parámetros . . . . . . . . . . . . . . . . . . . . . 17 2.5. Cadenas de Markov v́ıa Monte Carlo . . . . . . . . . . . . . . . . . . . 18 2.6. Estado actual del problema . . . . . . . . . . . . . . . . . . . . . . . . 19 3. Método de trabajo 22 3.1. Análisis de los datos simulados . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1. Datos simulados . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.2. Precisión de la estimación . . . . . . . . . . . . . . . . . . . . . 24 3.1.3. Capacidad para predecir . . . . . . . . . . . . . . . . . . . . . . 25 3.1.4. Tiempo computacional . . . . . . . . . . . . . . . . . . . . . . . 26 3.2. Análisis de los datos reales . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3. Implementación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.1. Software y Hardware . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.2. Obtención del número de part́ıculas . . . . . . . . . . . . . . . . 27 3.4. Algoritmos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4.1. Movimiento browniano fraccionario estándar . . . . . . . . . . . 28 3.4.2. Filtro de part́ıculas Liu-West . . . . . . . . . . . . . . . . . . . 29 3.4.3. Filtro de part́ıculas marginal Metropolis-Hastings . . . . . . . . 31 3.4.4. Filtro de part́ıculas secuencial aumentado MCMC . . . . . . . . 33 3.4.5. Filtro de part́ıculas secuencial aumentado MCMC con remuestreo 35 v 4. Resultados 37 4.1. Filtro de part́ıculas Liu-West . . . . . . . . . . . . . . . . . . . . . . . . 37 4.1.1. Sensibilidad al número de part́ıculas . . . . . . . . . . . . . . . 41 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings . . . . . . . . . . . . 42 4.2.1. Distribución a priori Gamma . . . . . . . . . . . . . . . . . . . 48 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC . . . . . . . . . . . 49 4.3.1. Sensibilidad al periodo de quema en los MCMC’s . . . . . . . . 55 4.3.2. Sensibilidad al periodo de quema y número de part́ıculas . . . . 56 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo . . . 58 4.5. Resumen de resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.6. Aplicación con datos reales . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.6.1. Datos de un mes . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.6.2. Datos de un año . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5. Conclusiones generales y recomendaciones 71 5.1. Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2. Recomendaciones de trabajos futuros . . . . . . . . . . . . . . . . . . . 74 Apéndices 75 A. Datos 76 A.1. Datos simulados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 A.2. Datos Reales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 B. Códigos 83 Referencias bibliográficas 99 vi RESUMEN En este estudio se trabaja con un modelo de volatilidad estocástica con memo- ria larga descrito por medio de un modelo de espacio estado. Se tienen dos procesos estocásticos, uno de ellos cuenta con información observada y el otro es un proceso latente. Las ecuaciones que describen estos procesos cuentan con 1 y 3 parámetros respectivamente. Estos parámetros y las variables del modelo se deben estimar. El problema se abordó aplicando una combinación del algoritmo de filtro de part́ıcu- las conocido como filtro de part́ıculas bootstrap y métodos de Cadenas de Marvok v́ıa Monte Carlo. Espećıficamente se aplicaron los algoritmos: filtro de part́ıculas Liu-West, filtro de part́ıculas marginal Metropolis-Hastings, filtro de part́ıculas secuencial aumen- tado MCMC y secuencial aumentado MCMC con remuestreo. Se obtuvo que el filtro de part́ıculas secuencial aumentado MCMC con remuestreo es el mejor algoritmo según las medidas usadas (Ráız del error cuadrático medio y el score de intervalo). Este algoritmo logra estimaciones buenas, tanto para los parámetros como para las variables del modelo. Además se ejecutó el algoritmo con datos reales usando el ı́ndice S&P500. vii LISTA DE CUADROS 3.1. Resultados de procedimiento para determinar el número de part́ıculas. . 28 4.1. Resumen de valores estad́ısticos de las muestras obtenidas con LW. . . 38 4.2. Valor promedio, desviación estándar y coeficiente de variación de la es- timación de los parámetros junto con el valor real a estimar con LW. . 38 4.3. Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 80 % y 90 % de las variables del modelo. Algoritmo LW. . . . . . . . . . 41 4.4. Coeficiente de variación de los escenarios de sensibilidad al número de part́ıculas en LW. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.5. Valores promedio y desviación estándar para diferentes escenarios de número de part́ıculas en LW. . . . . . . . . . . . . . . . . . . . . . . . . 42 4.6. Resumen de valores estad́ısticos de las muestras obtenidas con PMMH. 43 4.7. Valor promedio, desviación estándar y coeficiente de variación de la es- timación de los parámetros junto con el valor real a estimar con PMMH. 43 4.8. Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 80 % y 90 % de las variables del modelo. Algoritmo PMMH. . . . . . . 48 4.9. Resumen de valores estad́ısticos de las muestras obtenidas con PMMH Gamma. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 viii 4.10. Valor promedio, desviación estándar y coeficiente de variación de la es- timación de los parámetros junto con el valor real a estimar con PMMH Gamma. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.11. Resumen de valores estad́ısticos de las muestras obtenidas con SAMCMC. 51 4.12. Valor promedio y desviación estándar de la estimación de los parámetros junto con el valor real a estimar con SAMCMC. . . . . . . . . . . . . . 51 4.13. Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 90 % y 80 % de las variables del modelo. Algoritmo SAMCMC. . . . . . 54 4.14. Coeficiente de variación de los escenarios de sensibilidad al periodo de quema algoritmo SAMCMC. . . . . . . . . . . . . . . . . . . . . . . . . 55 4.15. Valores promedio y desviación estándar para diferentes escenarios de periodo de quema en los MCMC’s algoritmo SAMCMC. . . . . . . . . 56 4.16. Coeficiente de variación de los escenarios de sensibilidad al periodo de quema y número de part́ıculas. Algoritmo SAMCMC. . . . . . . . . . . 57 4.17. Valores promedio y desviación estándar para diferentes escenarios de periodo de quema y número de part́ıculas en los MCMC’s. Algoritmo SAMCMC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.18. Resumen de valores estad́ısticos de las muestras obtenidas con SAMCMC con remuestreo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.19. Valor promedio y desviación estándar de la estimación de los parámetros junto con el valor real a estimar con SAMCMC con remuestreo. . . . . 60 4.20. Ráız del error cuadrático medio (RECM) e Interval score (IS) con nive- les 90 % y 80 % de las variables del modelo. Algoritmo SAMCMC con remuestreo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.21. Comparación de la ráız del error cuadrático medio. Promedios sobre las 10 ejecuciones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.22. Comparación del score de intervalo. Promedios sobre las 10 ejecuciones. 65 ix 4.23. Comparación del score de intervalo. Promedios sobre las 10 ejecuciones. 65 4.24. Valores estad́ısticos de los parámetros obtenidos con SAMCMC con re- muestreo en datos de S&P500 del periodo 2008-12-29 al 2009-01-29. Los valores de Chronopoulou corresponden al estudio realizado en [6]. . . . 66 4.25. Valores estad́ısticos de los parámetros obtenidos con SAMCMC en datos de S&P500 del periodo 2008-12-29 al 2009-12-28. . . . . . . . . . . . . 70 A.1. Detalle de los valores simulados. . . . . . . . . . . . . . . . . . . . . . . 76 A.2. Información histórica del ı́ndice S&P500. . . . . . . . . . . . . . . . . . 79 x LISTA DE FIGURAS 3.1. Caminos simulados de las variables Log-Precio y volatilidad usados para ejecutar los algoritmos. . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.1. Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo LW . . . . . . . . . . . . . . 39 4.2. Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua ro- ja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Ĺıneas grises muestra de caminos estimados. Algoritmo LW . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.3. Histogramas de las muestras obtenidas con PMMH. La ĺınea discontinua de color rojo muestra el valor real. El área sombreada color gris corres- ponde a la estimación de la densidad por medio de un kernel gaussiano. 44 4.4. Gráfico de traza de las muestras obtenidas con PMMH. La ĺınea discon- tinua de color rojo muestra el valor real. . . . . . . . . . . . . . . . . . 45 4.5. Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo PMMH. . . . . . . . . . . . 46 xi 4.6. Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja valor promedio de la estimación, banda gris intervalo de confianza 90 %. Algoritmo PMMH. . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.7. Histogramas de las muestras obtenidas con SAMCMC. La ĺınea discon- tinua de color rojo muestra el valor real. El área sombreada color gris corresponde a la estimación de la densidad por medio de un kernel gaus- siano. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.8. Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo SAMCMC. . . . . . . . . . 52 4.9. Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja valor promedio de la estimación, banda gris intervalo de confian- za 90 %. Ĺıneas grises son muestras de caminos estimados. Algoritmo SAMCMC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.10. Histogramas de las muestras obtenidas con SAMCMC con remuestreo. La ĺınea discontinua de color rojo muestra el valor real. El área som- breada color gris corresponde a la estimación de la densidad por medio de un kernel gaussiano. . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.11. Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo SAMCMC con remuestreo. . 61 4.12. Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja valor promedio de la estimación, banda gris intervalo de confian- za 90 %. Ĺıneas grises son muestras de caminos estimados. Algoritmo SAMCMC con remuestreo. . . . . . . . . . . . . . . . . . . . . . . . . . 62 xii 4.13. Camino del Log S&P500 del periodo 2008-12-29 al 2009-01-29. Ĺınea discontinua roja valor promedio de la estimación. La ĺınea discontinua verde mediana, banda gris intervalo de confianza 90 %. . . . . . . . . . 67 4.14. Camino de la volatilidad estimada usando los datos de S&P500 del pe- riodo 2008-12-29 al 2009-01-29. Ĺınea discontinua roja valor promedio y la ĺınea discontinua verde mediana. . . . . . . . . . . . . . . . . . . . . 68 4.15. Camino del Log S&P500 del periodo 2008-12-29 al 2009-12-28. Ĺınea dis- continua roja valor promedio de la estimación. Ĺınea discontinua verde mediana, banda gris intervalo de confianza 90 %. . . . . . . . . . . . . . 69 4.16. Camino de la volatilidad estimada usando los datos de S&P500 del pe- riodo 2008-12-29 al 2009-12-28. Ĺınea rojo es el promedio y la ĺınea verde es la mediana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 xiii LISTA DE ALGORITMOS 2.1. Remuestreo secuencial por importancia o Filtro de part́ıculas . . . . . . 15 2.2. Filtro de part́ıculas bootstrap . . . . . . . . . . . . . . . . . . . . . . . 16 2.3. Filtro de part́ıculas auxiliar . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1. Filtro de part́ıculas Liu-West (LW) . . . . . . . . . . . . . . . . . . . . 31 3.2. Filtro de part́ıculas marginal Metropolis-Hastings (PMMH) . . . . . . . 32 3.3. Filtro de part́ıculas secuencial aumentado MCMC . . . . . . . . . . . . 34 3.4. Remuestreo en SAMCMC . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.5. Remuestreo sistemático . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 xiv 1 CAPÍTULO 1 INTRODUCCIÓN 1.1. Antecedentes y justificación En la matemática financiera se han dado grandes esfuerzos para desarrollar modelos que describan la dinámica de los precios de un activo financiero. Un avance importante fue realizado por Paul Samuelson en los años sesentas al modificar el modelo de Bache- lier y aśı desarrollar el proceso geométrico browniano para el precio de un activo [4], a saber, dPt = µPt dt+ σ Pt dBt. (1.1) En la fórmula, el valor µ es la tasa de retorno esperada, σ es la volatilidad del precio del activo y Bt es un movimiento Browniano estándar [18]. El modelo geométrico tiene ventajas sobre el modelo de Bachelier. La expresión (1.1) afirma que el retorno esperado y la volatilidad (o incertidumbre) tienen el mismo valor sin importar el precio inicial del activo, aspecto en lo cual falla en capturar el modelo de Bachelier [18]. Con el modelo de Bachelier los precios pueden tomar valores 1.1. Antecedentes y justificación 2 negativos con una probabilidad positiva [4]; además los incrementos en los precios siguen una distribución normal con esperanza µ y varianza σ2, pero hay evidencia con datos observados que las colas son mucho más pesadas que las de la distribución normal [4]. Como el proceso geométrico tiene como solución una función del tiempo exponencial perturbada aleatoriamente [25], supera los dos problemas mencionados. Además es considerado un modelo razonable, como una primera aproximación al precio real de un activo financiero [25]. Se ha observado que la volatilidad vaŕıa en el tiempo, por ejemplo se puede ver el caṕıtulo 19 de [18]. En el modelo de Black-scholes [18, 25], el cual asume que el precio del activo subyacente sigue la fórmula (1.1), se obtiene la llamada volatilidad impĺıcita la cual muestra estar relacionada con el precio de ejercicio del activo subyacente. Ésto es una contradicción a lo asumido por el modelo geométrico [6, 7, 18]. Por otra parte, estudios muestran que para algunos activos financieros existe una dependencia de su valor actual con valores pasados [6, 7]. Por los motivos anteriormente mencionados, se ha propuesto que la volatilidad sea descrita por un proceso estocástico y reemplazar la ecuación (1.1) por, dPt = µPt dt+ σ(Xt)Pt dBt (1.2) donde σ(x) es una función determińıstica y Xt es un proceso estocástico. En nuestro caso el proceso Xt será un proceso de Ornstein-Uhlenbeck fraccionario [5, 21]. Por lo tanto el modelo de volatilidad estocástica que se va a usar en este estudio es el siguiente: dPt = µPt dt+ σ(Xt)Pt dBt dXt = αXt dt+ β dBH t (1.3) El término BH t corresponde a un proceso estocástico browniano fraccionario, si el parámetro H ∈ ]0.5, 1] el proceso es de memoria larga. Esto significa que el valor del proceso estocástico está fuertemente correlacionado con valores pasados, es decir la 1.1. Antecedentes y justificación 3 dependencia en términos de autocorrerlación decae muy lentamente [6]. Por otro lado, se tiene que esta variable aleatoria Xt no tiene observaciones expĺıcitas, es decir es una variable latente. En la matemática aplicada gran parte de los esfuerzos se enfocan en la estimación de los parámetros de los modelos desarrollados. Buenas aproximaciones de los parámetros producen resultados más confiables y precisos de los modelos. Este proyecto se centra en la estimación conjunta de los parámetros del modelo de volatilidad estocástica con memoria larga (1.3), incluyendo H el parámetro del proceso browniano fraccionario. En la literatura se encuentran art́ıculos donde se ha trabajado con el modelo (1.3). En el trabajo de Chronopoulou y Viens [6] trataron el problema de obtener el precio de opciones, realizaron un análisis de sensibilidad en la influencia de H sobre los resultados finales. Mostraron un método para estimar el parámetro de memoria larga, de forma separada del resto de parámetros, esto bajo un método impĺıcito en el cual con los datos generaron filtros de part́ıculas para un rango de valores de H para luego escoger el mejor bajo un criterio de error cuadrático medio [6]. En el art́ıculo [7], segunda parte del anterior, Chronopoulou y Viens trabajan el modelo (1.3) con tres versiones, continua, discretización de la versión continua y una versión discreta, pero usan el mismo método para estimar H descrito en [6]. En este trabajo se usó filtro de part́ıculas para abordar el problema de la estima- ción conjunta de los parámetros de (1.3). Los filtros de part́ıculas son una forma de Monte Carlo secuencial, basado en la técnica de muestreo llamada en inglés sampling importance resampling [28]. En su versión original el filtro de part́ıculas, es usado en modelos de espacio-estado1 para estimar la trayectoria o evolución de la variable latente. Supone que los parámetros del modelo son conocidos, uno de estos algoritmos es el filtro Bootstrap (FB) [12]. Han habido modificaciones en estos algoritmos con el fin de estimar conjuntamente la 1El modelo (1.3) es un modelo de espacio-estado. La primera ecuación corresponde a la variable observada y la segunda a la variable latente o no observada. 1.2. Objetivos 4 variable latente y los parámetros, tal es el caso del algoritmo de Liu-West (LW) [24]. El algoritmo de LW introduce un kernel suavisador en la dinámica del muestreo de los parámetros. Hay otros algoritmos que incorporan los métodos de muestreo de Monte Carlo via cadenas de Markov (MCMC por sus siglas en inglés). Este estudio trabaja con dos de estos algoritmos [19] y [20]. 1.2. Objetivos Los objetivos generales y espećıficos de este trabajo final de graduación son los siguientes: 1.2.1. Objetivo general Estimar de forma conjunta los parámetros del modelo de volatilidad estocástica a través de algoritmos de filtros de part́ıculas y MCMC. 1.2.2. Objetivos espećıficos • Implementar el algoritmo de Liu-West para estimar de manera conjunta los parámetros del modelo (1.3). • Implementar algoritmos que combinen filtro de part́ıculas y MCMC. • Comparar el desempeño de los algoritmos, usando para ello datos simulados y aplicar el mejor algoritmo a una muestra de datos reales. 5 CAPÍTULO 2 MARCO TEÓRICO 2.1. Proceso estocástico Partimos de un concepto fundamental dentro de la teoŕıa de la probabilidad, las variables aleatorias. Una variable aleatoria es una función que toma un evento, como puede ser el resultado de un experimento, y le otorga un valor numérico. Esta fun- ción es la base para construir herramientas más complejas y sofisticadas que permiten desarrollar modelos matemáticos para describir o pronosticar un fenómeno. Cuando el tiempo está involucrado en el fenómeno en el que se está interesado se puede utilizar el concepto de proceso estocástico. Se define un proceso estocástico como: Definición 2.1 ([25]) Un proceso estocástico X es una colección de variables aleato- rias (Xt, t ∈ T ) = (Xt(ω), t ∈ T, ω ∈ Ω), definido en algún espacio Ω. 2.1. Proceso estocástico 6 Una propiedad que es deseable en un proceso, es que su comportamiento en distri- bución no cambien cuando pasa el tiempo. Matemáticamente esto se define como: Definición 2.2 ([25]) Un proceso X es estrictamente estacionario si (Xt1 , . . . , Xtn) d = (Xt1+h, . . . , Xtn+h), para todo t1, . . . , tn ∈ T , n ≥ 0 y h tal que t1 + h, . . . , tn + h ∈ T . Con la siguiente definición sólo interesa que el comportamiento en distribución de los incrementos del proceso se mantengan sobre el tiempo. Un proceso que tenga esta propiedad no necesariamente es estacionario [25]. Definición 2.3 ([25]) Un proceso X se dice tener incrementos estacionarios si Xt −Xs d = Xt+h −Xs+h, para todo t, s ∈ T y h tal que t+ h, s+ h ∈ T . Además, X se dice tener incrementos independientes si para cualquier escogencia de ti ∈ T tal que t1 < · · · < tn y n ≥ 1, Xt2 −Xt1 , . . . , Xtn −Xtn−1 , son independientes. Una clase importante de procesos estocásticos son los auto-similares. Intuitivamente auto-similaridad significa que, patrones apropiadamente escalados de una trayectoria en cualquier intervalo de tiempo pequeño o grande tiene una forma similar [25]. Definición 2.4 ([2]) Se dice que Xt es un proceso H-auto-similar, si para algún factor 2.1. Proceso estocástico 7 c, el proceso rescalado c−HXc t, es igual en distribución al proceso original Xt. Es decir, (c−HXc t1 , . . . , c −HXc tk) d = (Xt1 , . . . , Xtk). Esta definición es necesaria para los procesos que se usaron en este trabajo. 2.1.1. Proceso estocástico de memoria larga Existen fenómenos que muestran ciertas caracteŕısticas, como lo explica Beran en [2], hay periodos relativamente largos donde las observaciones toman valores altos, pero también hay periodos largos donde las observaciones toman valores bajos. Además, la varianza de la muestra de observaciones parece tener promedio que decae a cero a una tasa más lenta que n−1 y la correlación ρ̂(k) decae a cero a una tasa proporcional a k−α, con α ∈ ]0, 1[. Todo esto nos lleva a la siguiente definición: Definición 2.5 ([2]) Sea Xt un proceso estacionario. Si existe un α ∈ ]0, 1[ y una constante c > 0 tal que ĺım k→∞ ρ(k) c k−α = 1. Entonces, Xt es llamado un proceso estacionario con memoria larga. El parámetro α se puede transformar y cambiar por H = 1 − α/2, y el concepto de memoria larga se logra cuando H ∈ ]0.5, 1[. Si un proceso Xt es H-auto-similar con incrementos estacionarios, entonces el pro- ceso Yt = Xt −Xt−1 tiene su función de autocorrelación dada por [2], ρ(k) = 1 2 [(k + 1)2H − 2k2H + (k − 1)2H ]. (2.1) Usando teoŕıa de series numéricas y desarrollo de Taylor, se llega a que ρ(k) H(2H − 1) k2H−2 −−−→ k→∞ 1, 2.1. Proceso estocástico 8 y cuando 0.5 < H < 1 se tiene que la autocorrelación decae muy lentamente a cero, de forma que: ∞∑ k=−∞ ρ(k) =∞. En resumen, un proceso H-auto-similar con incrementos estacionarios, sus incre- mentos tienen memoria larga si H ∈ ]0.5, 1[. 2.1.2. Movimiento browniano fraccionario y ruido gaussiano fraccionario Uno de los procesos estocásticos más utilizados es el movimiento browniano estándar, el cual es H-auto-similar con H = 1/2. El movimiento browniano estándar es definido como: Definición 2.6 ([25]) Un proceso estocástico B es llamado movimiento browniano (estándar) o proceso de Wiener si: • B0 = 0, • Bt tiene incrementos independientes y estacionarios, • Para todo t > 0, Bt tiene una distribución normal N(0, t), • Bt tiene trayectoria continua. Las trayectorias de un movimiento brownianoBt(ω), con ω fijo, son extremadamente irregulares a pesar de que hay continuidad. La razón principal de que esto suceda es porque sus incrementos son independientes [25]. En general, el movimiento browniano fraccionario estándar se define como: Definición 2.7 ([2, 7, 33]) Un movimiento browniano fraccionario estándar BH t es un proceso gaussiano H-auto-similar con H ∈ ]0, 1[ e incrementos estacionarios. Su 2.2. Modelo de volatilidad estocástica 9 distribución es definida por su función de covarianza Cov(BH t , B H s ) = 1 2 [t2H − (t− s)2H + s2H ]. (2.2) Ademas, Xt = BH k+1 - BH k es llamado ruido gaussiano fraccionario. Note que, si H = 1/2, BH t es un movimiento browniano estándar [2, 33]. Por otro lado, si H ∈ ]0.5, 1[ entonces Xt tiene memoria larga [2]. 2.2. Modelo de volatilidad estocástica Un modelo de volatilidad estocástica es un proceso estocástico en el cual su varianza o volatilidad es aleatoria [18]. Este tipo de modelo fue desarrollado con el fin de superar la debilidad del modelo de Black-Scholes. Al usar el proceso geométrico browniano para describir la dinámica del precio del activo subyacente se asume que la varianza es constante, pero la fórmula de Black-Scholes demuestra una relación entre la volatilidad y el precio de ejercicio del activo subyacente [7]. Lo anterior se puede deducir de las fórmulas para el precio de un contrato call europeo (2.3) y para el precio de un contrato put europeo (2.4) de donde se obtiene la volatilidad impĺıcita [18]: ct = N(d1)St −N(d2)K exp−r(T−t) (2.3) pt = N(−d2)K exp−r(T−t)−N(−d1)St (2.4) donde, d1 = ln(St/K) + (r + σ2/2)(T − t) σ √ T − t (2.5) d2 = d1 − σ √ T − t (2.6) 2.2. Modelo de volatilidad estocástica 10 y N es la función de distribución de una normal estándar, (T − t) es el tiempo al vencimiento, St es el precio spot, K es el precio de ejercicio, r es la tasa libre de riesgo y σ es la volatilidad del retorno del activo subyacente. También hay evidencia documentada de la llamada persistencia de la volatilidad [9, 6, 33]. La persistencia de la volatilidad significa que la volatilidad de hoy nos dice algo no sólo sobre la volatilidad de mañana, sino también sobre la volatilidad de más d́ıas en el futuro. En este sentido, Comte [9] introduce el concepto de memoria larga en el modelo continuo de volatilidad estocástica de Hull and White. En este estudio se trabajó con el modelo propuesto por Comte [9]. En este modelo la volatilidad se describe a través de un proceso estocástico, con el fin de solucionar el inconveniente descrito anteriormente. dPt = µPt dt+ σ(Xt)Pt dBt (2.7) dXt = −αXt dt+ β dBH t (2.8) Pero en vez de usar directamente el precio se usa el logaritmo del precio Yt = log(Pt), con lo cual el modelo es el siguiente: dYt = ( µ− σ(Xt) 2 2 ) dt+ σ(Xt) dBt (2.9) dXt = −αXt dt+ β dBH t (2.10) La ecuación (2.9) describe la dinámica del precio de un activo financiero, es un proceso geométrico browniano con la volatilidad definida por σ(x) una función deter- mińıstica y Xt un proceso estocástico. La dinámica de Yt se describe con la ecuación (2.10), un proceso de Ornstein-Uhlenbeck fraccionario [5, 21]. En los art́ıculos de [6] y [7] se utiliza este modelo. 2.3. Modelo espacio estado 11 2.3. Modelo espacio estado Un modelo de espacio estado describe la dependencia probabiĺıstica entre dos va- riables aleatorias a través del tiempo y para ello involucra dos series de tiempo. Una de las series (yt) cuenta con observaciones y la otra serie (xt) no tiene observaciones, en otras palabras es una serie de tiempo de variables latentes. El modelo define la relación entre las series de la siguiente forma [22]: xt = Ft(xt−1, εt) (2.11) yt = Ht(xt, εt) (2.12) donde εt y εt son errores no necesariamente gaussianos, F y G son funciones no lineales. La primer ecuación (2.11) llamada ecuación de estado, relaciona el estado actual del sistema con el estado anterior tomando en cuenta un componente de error. La ecuación (2.12) se le llama ecuación de observación y relaciona las observaciones con el estado actual del modelo agregando un error [4]. En los modelos de espacio estado el principal problema es la estimación del vector de estados x0:t = {x1, . . . , xt} a partir del vector de observaciones y0:t = {y1, . . . , yt} [22]. Lo anterior se puede hacer en un sentido bayesiano, estimando la distribución posterior de todos los estados dada todas las observaciones, esto es p(x0:t | y1:t) = p(y1:t |x0:t) p(x0:t) p(y1:t) . Cada vez que se obtiene una nueva observación se tiene que recalcular la distribución y cuando t es grande el cálculo prácticamente no se puede realizar [12, 32]. En lugar de estimar la distribución completa se puede trabajar con la distribución marginal del estado latente y con un modelo probabilistico con propiedades de Markov [12, 23, 32], aśı el modelo de espacio estado descrito por las ecuaciones (2.11) y (2.12) se puede 2.3. Modelo espacio estado 12 reescribir como, Probabilidad inicial x0 ∼ p(x0) Probabilidad de transición xk ∼ p(xt |xt−1), t ≥ 1 Probabilidad de la observación dado el estado yt ∼ p(yt |xt), t ≥ 1. El modelo va a tener las siguientes propiedades. Propiedad 2.1 ([32]) El estado actual xt dado xt−1 es condicionalmente indepen- diente de los estados pasados antes del tiempo t− 1, esto es, p(xt |x1:t−1, y1:t−1) = p(xt |xt−1). Propiedad 2.2 ([32]) La observación actual yt dado xt es condicionalmente indepen- diente de los estados y observaciones anteriores, esto es, p(yt |x1:t, y1:t−1) = p(yt |xt). Las distribuciones que se van a considerar ahora son [22]: Tipo Fórmula Filtro p(xt | y1:t) Predicción p(xt+n | y1:t) y usando la fórmula de Bayes y la ecuación de Chapman-Kolmogorov se obtienen las siguientes fórmulas [12, 23] para obtener las distribuciones, p(xt | y1:t−1) = ∫ p(xt|xt−1) p(xt−1 | y1:t−1) dxt−1 (2.13) p(xt | y1:t) = p(yt |xt) p(xt | y1:t−1)∫ p(yt |xt) p(xt | y1:t−1) dxt (2.14) 2.4. Filtro de part́ıculas 13 2.4. Filtro de part́ıculas Una de las formas de obtener una estimación de las distribuciones (2.13) (2.14) es por medio de simulación de Monte Carlo. La forma directa de aplicar Monte Carlo es obteniendo muestras de las distribuciones involucradas [17, 29, 30], en nuestro caso p(x | y1:t), pero no siempre es posible. En términos generales se quiere ser capaz de estimar E[g(x) | y1:t] = ∫ g(x) p(x | y1:t) dx. Otra alternativa para realizar la estimación es usando el método Muestreo por im- portancia [30, 31]. Se incluye una distribución llamada distribución de muestreo por importancia π(x | y1:t) de la cual se puede obtener muestras fácilmente [12]. El método está basado en la siguiente descomposición, en donde se asume que el soporte de π contiene el soporte de p(x | y1:t) [30], E[g(x) | y1:t] = ∫ g(x) p(x | y1:t) dx (2.15) = ∫ g(x) [ p(x | y1:t) π(x | y1:t) ] π(x | y1:t) dx (2.16) = ∫ g(x)ω(x) π(x | y1:t) dx. (2.17) Para aproximar (2.17) se utiliza simulación Monte Carlo, al obtener muestras de la distribución de muestreo por importancia: x(i) ∼ π(x | y1:t), i = 1, . . . , N 2.4. Filtro de part́ıculas 14 y la aproximación sigue de la forma [12, 32]: E[g(x) | y1:t] ≈ 1 N N∑ i=1 p(x(i) | y1:t) π(x(i) | y1:t) g(x(i)) (2.18) = N∑ i=1 ω̃(x(i)) g(x(i)). (2.19) La densidad posterior se estima con este método como [12], p(x | y1:t) ≈ N∑ i=1 ω̃(x(i)) δx(x (i)), (2.20) donde δ(x) es la función delta de Dirac. Como se menciona en [12] este método no es adecuado cuando se tiene que calcular la evolución secuencial de distribuciones posteriores, cada vez que una observación es disponible se tiene que recalcular los pesos y con el paso del tiempo se complican más los cálculos. Para poder sobrellevar el inconveniente anteriormente mencionado se desarrolla una versión llamada Muestreo secuencial por importancia [12, 22, 30, 32]. La modificación parte de la distribución completa posterior y de las propiedades 2.1 y 2.2 para tener una fórmula recursiva [12]: p(x0:t | y1:t) = p(x0:t−1 | y1:t−1) p(yt |xt) p(xt |xt−1) p(yt | y1:t−1) (2.21) ∝ p(x0:t−1 | y1:t−1) p(yt |xt) p(xt |xt−1). (2.22) Si se escoge la distribución de muestreo por importancia tal que cumpla: π(x0:t | y1:t) = π(xt |x0:t−1, y1:t) π(x0:t−1 | y1:t−1), 2.4. Filtro de part́ıculas 15 se tiene que los pesos también se pueden calcular de forma recursiva, ω (i) t ∝ p(x (i) 0:t−1 | y1:t−1) π(x (i) 0:t−1 | y1:t−1) p(yt |x(i) t ) p(x (i) t |x (i) t−1) π(x (i) t |x (i) 0:t−1, y1:t) (2.23) = ω (i) t−1 p(yt |x(i) t ) p(x (i) t |x (i) t−1) π(x (i) t |x (i) 0:t−1, y1:t) . (2.24) De esta forma sólo hay que obtener muestras de xt e ir construyendo x0:t = {x0:t−1, xt}. Algoritmo 2.1 Remuestreo secuencial por importancia o Filtro de part́ıculas 1: Obtener muestras x (i) 0 de la distribución a priori y tomar ω (i) 0 = 1/N , x (i) 0 ∼ p(x0), i = 1, . . . , N. 2: Para t = 1, . . . , T hacer 3: Obtener muestras x (i) t de la distribución por importancia x (i) t ∼ π(xt |x(i) t−1, y1:t), i = 1, . . . , N. 4: Calcular los pesos nuevos ω (i) t ∝ p(yt |x(i) t ) p(x (i) t |x (i) t−1) π(x (i) t |x (i) t−1, y1:t) , y normalizarlos para que sumen 1. 5: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t : i = 1, . . . , N} con probabilidades {ω(i) t : i = 1, . . . , N}, luego tomar ω (i) t = 1/N . 6: Fin Para El muestreo secuencial por importancia también tiene sus debilidades, hay que ir almacenando la historia de los estados x0:t para poder aplicar el algoritmo [32]. Esto se resuelve escogiendo la distribución de muestreo por importancia de forma conveniente [32], π(xt |x0:t−1, y1:t) = π(xt |xt−1, y1:t). 2.4. Filtro de part́ıculas 16 Además se da el problema de degeneración [30], fácilmente se puede encontrar la si- tuación de que después de un tiempo casi todas las part́ıculas tienen pesos con valor muy cercano a cero. Para solucionar el problema de la degeneración se incorpora un remuestreo con el fin de eliminar las part́ıculas con pesos muy pequeños y duplicar las part́ıculas con pesos grandes, esto da lugar a la versión de muestreo por importan- cia llamado Remuestreo secuencial por importancia también conocido como Filtro de part́ıculas [12, 22]. El Algoritmo 2.1 resume el método. Algoritmo 2.2 Filtro de part́ıculas bootstrap 1: Obtener muestras x (i) 0 de la distribución a priori y tomar ω (i) 0 = 1/N , x (i) 0 ∼ p(x0), i = 1, . . . , N. 2: Para t = 1, . . . , T hacer 3: Obtener muestras x (i) t de la distribución por importancia x (i) t ∼ p(xt |x(i) t−1), i = 1, . . . , N. 4: Calcular los pesos nuevos ω (i) t ∝ p(yt |x(i) t ), y normalizarlos para que sumen 1. 5: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t : i = 1, . . . , N} con probabilidades {ω(i) t : i = 1, . . . , N}, luego tomar ω (i) t = 1/N . 6: Fin Para Dentro de las variaciones al algoritmo de filtro de part́ıculas la más conocida es el Filtro de part́ıculas bootstrap. Usa como distribución de muestreo por importancia la distribución de transición p(xt |xt−1) del modelo, con esto los pesos se obtienen fácil- mente, ω (i) t ∝ p(yt |x(i) t ). La implementación puede ser muy fácil pero puede requerir un gran número de part́ıculas para tener buenos resultados [32]. El Algoritmo 2.2 resume el método. Otra variación importante, es el llamado Filtro de part́ıculas auxiliar [26]. Realiza un 2.4. Filtro de part́ıculas 17 remuestreo de las part́ıculas del paso anterior con pesos proporcionales a p(yt | g(x (i) t−1)), donde g puede ser la esperanza o moda de p(xt |xt−1). La idea es tener las part́ıculas con mayor probabilidad de que su propagación esté en la zona de la densidad posterior. El Algoritmo 2.3 resume el método. Algoritmo 2.3 Filtro de part́ıculas auxiliar 1: Obtener muestras x (i) 0 ∼ p(x0) de la distribución a priori y tomar ω (i) 0 = 1/N, i = 1, . . . , N. 2: Para t = 1, . . . , T hacer 3: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t−1 : i = 1, . . . , N} con probabilidades {w(i) t = p(yt | g(x (i) t−1)) : i = 1, . . . , N}. 4: Obtener muestras x (i) t de la distribución por importancia, x (i) t ∼ p(xt | x̃(i) t−1), i = 1, . . . , N. 5: Calcular los pesos nuevos ω (i) t ∝ p(yt |x(i) t ) p(yt | g(x̃ (i) t−1)) , y normalizarlos para que sumen 1. 6: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t : i = 1, . . . , N} con probabilidades {ω(i) t : i = 1, . . . , N}, luego tomar ω (i) t = 1/N . 7: Fin Para 2.4.1. Estimación de parámetros El modelo de espacio estado (2.11), (2.12) por lo general tiene parámetros θ no conocidos que se deben estimar, aśı incorporar los parámetros permite que el modelo 2.5. Cadenas de Markov v́ıa Monte Carlo 18 pueda ser ajustado a datos reales [32]: Probabilidad inicial de los parámetros θ ∼ p(θ) Probabilidad inicial del sistema x0 ∼ p(x0 |θ) Probabilidad de transición xk ∼ p(xt |xt−1,θ), t ≥ 1 Probabilidad de la observación dado el estado yt ∼ p(yt |xt,θ), t ≥ 1. Y la distribución que estamos interesados en estimar es: p(xt,θ | y1:t) ∝ p(yt |xt,θ) p(xt | y1:t−1,θ) p(θ | y1:t−1). 2.5. Cadenas de Markov v́ıa Monte Carlo Las cadenas de Markov v́ıa Monter Carlo (MCMC por sus siglas en inglés) son métodos que generan muestras para θ provenientes de la distribución posterior p(θ | y) cuando esta distribución no puede ser directamente simulada [29]. Los métodos MCMC son Monte Carlo porque dependen de muestreos aleatorios. Producen cadenas de Markov pues las muestras que se generan no son independientes, cada muestra depende de la anterior. Además su distribución estacionaria es p(θ | y) [30]. Un mecanismo para la generación de una cadena de este tipo es llamado algoritmo de Metropolis-Hastings el cual se puede describir de la siguiente manera: 1. Generar θ∗ ∼ J(θ |θ(s)) 2. Calcular la tasa de aceptación α = p(θ∗ | y) p(θ(s) | y) = p(y |θ∗) p(θ∗) p(y |θ(s)) p(θ(s)) 2.6. Estado actual del problema 19 3. Hacer θ(s+1) = θ ∗ con probabilidad mı́n(α, 1) θ(s) con probabilidad 1−mı́n(α, 1). La distribución p(θ | y) es llamada la distribución objetivo y J(θ |θ(s)) es llamada la distribución propuesta. La propuesta más común de usar (desde que Hastings propuso su uso [30]) es la caminata aleatoria θ∗ = θ(s) + ε y las distribuciones más usadas para la perturbación ε suelen ser uniforme, normal o Cauchy [29]. En la práctica se suele no considerar las primeras B iteraciones del proceso, es decir sólo se considera la muestra {θB+1, . . . ,θm}. Estas B iteraciones iniciales se le llama periodo de quema y el propósito es no considerar el periodo en que la cadena se mueve desde su valor inicial hasta la zona del espacio de parámetros con mayor probabilidad posterior [17]. 2.6. Estado actual del problema El problema que se trata en este estudio tiene sus antecedentes en los art́ıculos de Chronopoulou [6, 7], en el cual el problema de interés fue estimar el precio de opciones considerando el modelo (2.9), (2.10) para el precio de un activo subyacente. Parte del art́ıculo [6] describió como estimar los parámetros del modelo incluyendo el parámetro de memoria larga. La estimación no se realizó de forma conjunta, se usó un método para estimar H para luego estimar los otros parámetros del modelo (2.9), (2.10). La forma de obtener H es la siguiente: toman valores de H desde 0.5 hasta 0.95 con incrementos de 0.01, para cada uno de los valores estiman la distribución de la volatilidad, calculan el 2.6. Estado actual del problema 20 correspondiente precio de la opción para varios precios de ejercicio, para luego obtener el error cuadrático medio del precio de la opción con respecto al bid-ask spread. Se escoge el H correspondiente al valor más pequeño del error cuadrático medio. Con respecto al método usado en este trabajo, filtro de part́ıculas, se tiene conoci- miento de varios art́ıculos con información valiosa. Recientemente se ha incorporado en los filtros de part́ıculas las cópulas con el fin de estimar los parámetros de un modelo de espacio estado. El art́ıculo de Deng [10] desarrolla un modelo de espacio estado en el cual la varianza del error es función del estado del sistema, además modifica el filtro Liu-West incorporandole cópulas para considerar la dependencia de los parámetros con la variable del sistema y las observaciones. Se cambia la forma de calcular los pesos, se usa la función de densidad de la cópula para calcularlos: w (i) t = c(yt | g(x (i) t−1),m (i) t−1) ω (i) t ∝ c(yt |x(i) t ,θ (i) t ) c(yt | g(x̃ (i) t−1), m̃ (i) t−1) , donde c es la función de densidad de la cópula. El art́ıculo concluye, basado en los resultados numéricos del error cuadrático medio, que el filtro de part́ıculas propuesto es superior al filtro de Liu-West. El art́ıculo de Fan [13] desarrolla un complicado filtro de part́ıculas incluyendo cópulas, el cual lo aplica a problemas de hidroloǵıa. El trabajo con cópulas, además de ser usadas para incluir las estructuras de dependencia, se hace para corregir dos problemas que se dan con los filtros de part́ıculas, a saber, la degeneración de las part́ıculas y el empobrecimiento de las muestras. Se concluye que el método propuesto tiene más exactitud que los filtros tradicionales usando muestras de tamaño de 100 part́ıculas y que el tamaño de la muestra de part́ıculas no afecta significativamente el desempeño del algoritmo propuesto. Un art́ıculo importante es el de Jacob [19], en donde se trabajó en el campo de 2.6. Estado actual del problema 21 la electroqúımica. Usan un modelo de espacio estado donde la variable latente no es Markov. Los autores indican que los métodos de filtros y filtros MCMC se pueden implementar directamente en procesos que no cumplen la propiedad de Markov. El problema de la estimación de parámetros se aborda con un método de cadenas de Markov v́ıa Monte Carlo1 llamado Filtro de part́ıculas marginal Metropolis-Hastings. El art́ıculo estudia la identificación de los parámetros y su sensibilidad a la escogencia de distribuciones a priori y al número de observaciones. Con respecto a la inferencia bayesiana se sugiere probar con otros métodos como Filtro de part́ıculas Gibbs y SMC2. Otro art́ıculo importante para este estudio es el de Javvad [20]. En el art́ıculo se implementa un algoritmo denominado secuencial aumentado de cadenas de Markov v́ıa Monte Carlo para obtener la distribución posterior de los parámetros. El algoritmo se desarrolla mejorando los métodos secuenciales actuales en los que se debe contar con estad́ısticos suficientes o el uso de un modelo artificial para la evolución de los parámetros. 1MCMC por sus siglas en inglés. 22 CAPÍTULO 3 MÉTODO DE TRABAJO El estudio se llevará a cabo con dos tipos de datos, datos simulados y datos reales. Los datos simulados nos permiten identificar debilidades y beneficios de cada uno de los algoritmos, aśı como realizar comparaciones entre los diferentes métodos, mientras que con los datos reales se puede confirmar la aplicabilidad de los filtros de part́ıculas para un problema en espećıfico [28]. Los algoritmos se implementaron en el lenguaje estad́ıstico R [27]. Por otro lado, se usa el promedio como estimador de los parámetros y variables objetivos del estudio. Este caṕıtulo es organizado de la siguiente forma: en la sección 3.1 se indica como se obtuvieron los datos simulados, además de las medidas usadas para medir el des- empeño de los algoritmos. La sección 3.2 indica cuales datos reales fueron usados. La implementación de los algoritmos se trata en la sección 3.3, y finalmente en la sección 3.4 se detallan los algoritmos que se usaron para resolver el problema. 3.1. Análisis de los datos simulados 23 3.1. Análisis de los datos simulados 3.1.1. Datos simulados Para medir el rendimiento de los algoritmos se utilizaron datos simulados de medi- ciones diarias en un periodo de tiempo, utilizando el modelo (2.9), (2.10). Para esto se usó la discretización de Euler de primer orden [7], a saber, Yt+1 = Yt + (µ−X2 t /2) dt+Xt εt √ dt (3.1) Xt+1 = Xt − αXt dt+ β (BH t+1 −BH t ). (3.2) Donde, Yt = log(Pt) es el logaritmo del precio del stock, Xt es la volatilidad, εt es el error que sigue una distribución normal estándar, dt = 1/N y BH t es un movimiento browniano fraccionario estándar. Se tomaron valores espećıficos para los parámetros, a saber: α = 0.02733, β = 0.07567, µ = 0.0014, H = 0.6, x0 = 6.802, y0 = 0.35 y se tomó σ(x) = x, función que se usa en la ecuación (2.9). A excepción de µ, estos valores se usaron en [6]. Además se definieron valores para N (el número de part́ıculas) y T (cantidad de datos históricos), de 130 y 255 respectivamente. Las series simuladas de las variables se pueden observar en la Figura 3.1. En términos generales se realizaron pruebas de sensibilidad en los algoritmos, con el fin de observar el efecto en los resultados. En el algoritmo LW se varió el número de part́ıculas. En PMMH se cambió la distribución a priori con el objetivo de ver la influencia en los resultados, se varió de una distribución normal multivariada a una gamma multivariada. Para el algoritmo SAMCMC se varió el periodo de quema en los MCMC’s y el número de part́ıculas. Por otro lado, se realizaron mediciones del desempeño tales como: precisión, capacidad para predecir y tiempo de ejecución de los diferentes algoritmos. Con estas medidas se determinó cual de los algoritmos es mejor. 3.1. Análisis de los datos simulados 24 Los algoritmos se ejecutaron diez veces (K = 10), con el objetivo de aislar el efecto de la aleatoriedad de la simulación en las medidas de desempeño. 0.34 0.36 0.38 0.40 0 100 200 Día V ol at ili da d Volatilidad subyacente simulada (H = 0.6) 6.6 6.7 6.8 6.9 7.0 7.1 0 100 200 Día P re ci o S to ck Log−Precios de stock simulado Figura 3.1: Caminos simulados de las variables Log-Precio y volatilidad usados para ejecutar los algoritmos. 3.1.2. Precisión de la estimación Para determinar la precisión de la estimación, esto es medir la desviación entre el valor real y el valor estimado, se calculó para cada K la ráız del error cuadrático medio (RECM) y el valor que se reporta es le promedio sobre K [13, 20, 28]. La medida se calcula para las variables del modelo espacio estado Xt y Yt, y también 3.1. Análisis de los datos simulados 25 para los parámetros del modelo α, β, H y µ. Recordemos que lo que se busca en la simulación es analizar la capacidad predictiva de los estimadores en parámetros y en la variable latente. RECM = 1 K K∑ k=1 RECMk = 1 K K∑ k=1 1 T T∑ t=1 √√√√ N∑ j=1 (xt − X̂j,t,k) 2 N . (3.3) RECM = 1 K K∑ k=1 RECMk = 1 K K∑ k=1 √√√√ N∑ j=1 (θ − θ̂j,k)2 N . (3.4) Donde, xt es el valor observado simulado de las variables del modelo espacio estado X̂j,t,k es la muestra j del valor estimado para el momento t en la ejecución k y T es largo del camino observado. θ es el valor del parámetro, θ̂j,k es el valor estimado y N es el tamaño de la muestra. Un valor de RECM cercano a cero refleja en forma total la cercańıa del estimador al valor real e indica que el algoritmo en estudio es efectivo en inferir la variable en el modelo [20]. 3.1.3. Capacidad para predecir El score de intervalo [14] es un indicador diseñado para predicciones en formato de intervalos de percentiles, este indicador sólo requiere del intervalo de predicción central (1− α)× 100 %. La definición del score es: IS1−α = (u− l) + 2 α (l − x)1{xu}. (3.5) Donde 1 es la función indicadora, l es el percentil α/2 y u es el percentil 1 − α/2 de la distribución emṕırica posterior de la variable del modelo espacio estado o del parámetro. El objetivo es medir la capacidad predictiva del estimador, al recompensar los inter- 3.2. Análisis de los datos reales 26 valos de predicción más angostos y penalizar si el valor real queda afuera del intervalo de predicción central. Note que entre más ancho el intervalo de predicción central más severa es la penalización [3]. Se reporta para la estimación de los parámetros el promedio de los score de intervalo de las K ejecuciones. IS1−α = 1 K K∑ k=1 IS1−α,k (3.6) Para la estimación de las variables del modelo de espacio estado se usa la siguiente fórmula: IS1−α = 1 K K∑ k=1 T∑ t=1 IS1−α,t,k T (3.7) donde T es el tamaño del del camino. Se calculó la medida con valores de α iguales a 0.1 y 0.2. 3.1.4. Tiempo computacional El tiempo de ejecución de los algoritmos se midió en segundos. Se reporta el pro- medio de las K ejecuciones. Si la duración es muy larga el resultado se muestra en horas. 3.2. Análisis de los datos reales Para el segundo tipo de datos se usó información real, espećıficamente el ı́ndice Standard & Poor’s 500 (S&P500) [34], uno de los más importantes de Estados Unidos. El ı́ndice será usado en la ejecución del algoritmo con mejor desempeño. Se tomó una serie corta y una larga del ı́ndice. Para la serie corta se tomaron datos entre las fechas 2008-12-29 al 2009-01-29, para un total de 22 observaciones. Para la serie larga se tomaron 252 observaciones del periodo de 2008-12-29 al 2009-12-28. Los datos se tomaron de Yahoo Finanzas y se presentan en el anexo A.2. 3.3. Implementación 27 3.3. Implementación En esta sección se detallan los temas involucrados en la implementación de los métodos. 3.3.1. Software y Hardware Los algoritmos se ejecutaron en una computadora con las siguientes caracteŕısticas: arquitectura x86 64, procesador Inter (R) Xeon (R) CPU E5-2630 v3 2.40GHZ, 16 núcleos, 62 Gb memoria RAM. Se usó el lenguaje estad́ıstico R [27] para programar y ejecutar los algoritmos en un sistema Linux Ubuntu 18.04.5 LTS. La versión de R que se utilizó fue 3.6.3. 3.3.2. Obtención del número de part́ıculas Para determinar un número de part́ıculas se siguió el procedimiento descrito en [19]. La estrategia consiste en tomar una muestra de la distribución a priori de θ, se ejecuta el filtro de part́ıculas 100 veces con el mismo θ. Se genera una serie de valores Z1, . . . , Z100 de p(y0:T |θ). Luego, se imita un Metropolis-Hastings usando como probabilidad α = mı́n(1, Zj/Zactual) para aceptar Zactual = Zj. Al final, se obtiene la tasa de veces que un nuevo Z fue aceptado. Según [19] una tasa de 10 % se puede considerar buena. Este procedimiento se ejecutó con varios valores de θ0 y número de part́ıculas. Los resultados se pueden ver en el cuadro 3.1. Se tomó la decisión de usar 130 part́ıculas como escenario base debido a que con las pruebas que se realizaron con 130 se tiene más de 10 % de tasa de aceptación. 3.4. Algoritmos 28 Cuadro 3.1: Resultados de procedimiento para determinar el número de part́ıculas. θ0 = (α0, β0, H0, µ0) Part́ıculas Tasa (0.06715, 0.05978, 0.78726, 0.00768) 50 15 % (0.05298, 0.06798, 0.26147, 0.00467) 130 11 % (0.01295, 0.07005, 0.41802, 0.00505) 100 11 % (0.02344, 0.05909, 0.60494, 0.00700) 100 10 % (0.05361, 0.04691, 0.97870, 0.00157) 130 20 % 3.4. Algoritmos En esta sección se describe los algoritmos que se implementaron, en particular: filtro de part́ıculas Liu-West, filtro de part́ıculas marginal Metropolis-Hastings, filtro de part́ıculas secuencial aumentado MCMC y filtro de part́ıculas secuencial aumentado MCMC con remuestreo. Además de estos algoritmos, se indica el algoritmo usado para la generación de los caminos del movimiento browniano fraccionario estándar. 3.4.1. Movimiento browniano fraccionario estándar Para generar caminos del movimiento browniano fraccionario estándar BH , se usa la descomposición de Choleski para la matriz de covarianza (2.2), (Γ)i,j = cov ( i N , j N ) , para i, j = 0, . . . , N − 1. (3.8) Donde Γ es simétrica definida positiva y admite la descomposición de Choleski, Γ = LLt con L matriz triangular inferior. Para simular un camino de BH , se genera un vector Z de una distribución normal estándar multivariada, se calcula LZ y la estimación del camino es B̃H = (0, (LZ)t)t 3.4. Algoritmos 29 [8] (al igual que en la definición 2.6 se asume que el proceso inicia en 0). Este método es exacto, pero tiene una complejidad computacional de O(N3) [8]. 3.4.2. Filtro de part́ıculas Liu-West Una técnica para poder incorporar la estimación de parámetros en los Algoritmos 2.2 y 2.3 es la llamada Kernel suavisador (KS), en donde se asume que la distribución de θ es una mezcla de normales multivariadas [24]. Si se tiene un conjunto de part́ıculas {x(i) t−1, θ (i) t−1 : i = 1, . . . , N} con sus respecti- vos pesos {ω(i) t−1 : i = 1, . . . , N} que aproximan a p(xt−1,θ | y1:t−1), se tendŕıa que la distribución posterior de θ se puede aproximar de la siguiente forma: p(θ | y1:t−1) ≈ N∑ i=1 ω (i) t−1N(θ |m(i) t−1, h 2Vt−1), donde, m (i) t−1 = aθ (i) t−1 + (1− a)θ, θ = N∑ i=1 θ (i) t−1 N , Vt−1 = N∑ i=1 [θ (i) t−1 − θ][θ (i) t−1 − θ]′ N , N(θ |m,V ) es una distribución normal multivariada, con esperanza m y matriz de covarianzas V . Además la constante a2 = 1− h2 es un factor de contracción [24]. Combinando el Filtro de part́ıculas auxiliar con kernel suavisador se obtiene el Filtro de Liu and West (LW) [23, 24, 28]. El Algoritmo 3.1 resume el método. Para aplicar este algoritmo se usó un factor de contracción, recomendado en [24] de a = 0.99495. La función g(x) que se usó es la esperanza. Los valores iniciales x0, y0 3.4. Algoritmos 30 y θ0 se obtienen de distribuciones normales multivariadas: x0 ∼ N(µx, σx,mı́n = 0,máx = 1), (3.9) donde, µx = sd(∇yobs)/ √ 1/N , y σx = 1. y0 ∼ N(µy, σy,mı́n = 0,máx =∞) (3.10) donde, µy = yobs, σy = var(yobs). θ0 ∼ N(µθ,Σθ,mı́n = (0, 0, 0, 0),máx = (1, 1, 1, 1)), (3.11) donde, µθ = (0, 0, 0.75, 0), Σθ = diag(5 · 10−4, 1 · 10−3, 5 · 10−2, 15 · 10−7). La escogencia de estos valores no afecta la calidad del ajuste del filtro de Liu-West. Para obtener los pesos w se utilizó la medida de probabilidad discreta usada en [6], la cual tiene la fórmula: p(yt |x(i) t ,θ (i) t ) = 3 √ n exp(−2|(yit − yt) 3 √ n|) n∑ i=1 3 √ n exp(−2|(yit − yt) 3 √ n|) , (3.12) donde, yit es la part́ıcula de la variable observada, yt es el valor observado y n es el número de part́ıculas. 3.4. Algoritmos 31 Algoritmo 3.1 Filtro de part́ıculas Liu-West (LW) 1: Para t = 1, . . . , T hacer 2: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t−1,θ (i) t−1 : i = 1, . . . , N} con probabilidades {w(i) t = p(yt | g(x (i) t−1),m (i) t−1) : i = 1, . . . , N}. 3: Obtener muestras θ (i) t del KS, θ (i) t ∼ N(θ | m̃(i) t−1, h 2Vt−1) 4: Obtener muestras x (i) t de la distribución por importancia, x (i) t ∼ p(xt | x̃(i) t−1,θ (i) t ), i = 1, . . . , N. 5: Calcular los pesos nuevos ω (i) t ∝ p(yt |x(i) t ,θ (i) t ) p(yt | g(x̃ (i) t−1), m̃ (i) t−1) , y normalizarlos para que sumen 1. 6: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t ,θ (i) t : i = 1, . . . , N} con probabilidades {ω(i) t : i = 1, . . . , N}, luego tomar ω (i) t = 1/N . 7: Fin Para 3.4.3. Filtro de part́ıculas marginal Metropolis-Hastings Uno de los algoritmos que se consideró fue el filtro de part́ıculas marginal Metropolis- Hastings [1]. Este algoritmo aproxima la distribución posterior p(θ | y0:T ) usando la combinación de los métodos: el filtro de part́ıculas [12] y las cadenas de Markov v́ıa Monte Carlo [17, 30]. El primero, es empleado para aproximar p(xT | y0:T ,θ) y la vero- similitud p(y0:T |θ). El segundo, se usa para aproximar p(θ | y0:T ) basado en la verosi- militud obtenida por el filtro de part́ıculas [11]. 3.4. Algoritmos 32 Algoritmo 3.2 Filtro de part́ıculas marginal Metropolis-Hastings (PMMH) 1: Obtener muestra de la distribución apriori: θ0 ∼ p(θ). 2: Ejecutar el filtro de part́ıculas, Algoritmo 2.2, con N part́ıculas para obtener p̃(y0:T |θ0). 3: Para m = 1, . . . ,M hacer 4: Propuesta: θ∗ ∼ q(θ |θm−1). 5: Ejecutar el filtro de part́ıculas, Algoritmo 2.2, con N part́ıculas para obtener p̃(y0:T |θ∗). 6: Calcular: α = mı́n ( 1, p(θ∗) p̃(y0:T |θ∗) p(θm−1) p̃(y0:T |θm−1) ) . 7: Hacer θm = { θ∗ con probabilidad mı́n(α, 1) θt−1 con probabilidad 1−mı́n(α, 1). 8: Fin Para La expresión para la tasa α (tasa de aceptación), en el Algoritmo 3.2 punto 6, sugiere que el algoritmo efectivamente aproxima la densidad marginal p(θ | y0:T ), justificando el termino marginal [1]. El atractivo de este algoritmo es que la dificultad de muestrear de p(x0:T , θ | y0:T ) es reducida a muestrear en p(θ | y0:T ) [1]. Queda claro que este algoritmo se enfoca en estimar los parámetros del modelo en vez de estimar la variable latente. Para aplicar el filtro bootstrap, Algoritmo 2.2, se usó las distribuciones (3.9) y (3.10) para los valores iniciales. Para obtener los pesos w se utilizó la medida de probabilidad normal con promedio yt y desviación estándar 1, esto es w (i) t ∼ N(yt, 1). Se usó esta probabilidad porque se obtuvieron mejores resultados que al usar la pro- 3.4. Algoritmos 33 babilidad dada por (3.12). Para la distribución a priori y la propuesta se utilizaron distribuciones normales multivariadas truncadas, con valor mı́nimo (0, 0, 0, 0) y valor máximo (1, 1, 1, 1). La distribución a priori usa la esperanza (0, 0, 0.75, 0) y la matriz de covarianza es una matriz diagonal con elementos (0.001, 0.01, 0.01, 0.0000015). Para la distribución pro- puesta se usa una matriz de covarianza, una matriz diagonal con valores (0.00001, 0.001, 0.001, 0.0000001). Se decidió usar estos valores pues al usar otros valo- res los resultados del ajuste no vaŕıan mucho. 3.4.4. Filtro de part́ıculas secuencial aumentado MCMC Se implementó el algoritmo propuesto en [20]. El algoritmo es una combinación de los métodos: • Secuenciales, como LW, en donde el muestreo permite actualizaciones sucesivas de la distribución posterior a medida que hay nuevas observaciones disponibles, • y los métodos que usan MCMC, que permiten obtener muestras de la distribución objetivo posterior de θ, además que permiten obtener propuestas arbitrarias sin tener que asumir la existencia de estad́ısticos suficientes o un modelo artificial para la evolución de θ. El algoritmo propuesto en [20] construye muestras por medio de MCMC para cada t = 1, . . . , T , en las cuales sólo incorpora el t-ésimo componente p(yt | y1:t−1,θ) de la verosimilitud. Además, asume que la distribución posterior p(θ | y0:t) puede ser apro- ximada por una mixtura gaussiana. La mixtura gaussiana es obtenida a partir de la muestra (1θ B t , . . . , Mθ B t ) recolectada después del periodo de quema. El art́ıculo [20] realiza el ajuste para (θBt , xt). Para determinar el número de componentes de la mixtura se realizaron pruebas con diferentes valores. Se obtuvieron resultados buenos con un componente, es decir con 3.4. Algoritmos 34 una normal multivariada truncada con parámetros promedio y la matriz de covarianza de la muestra. Algoritmo 3.3 Filtro de part́ıculas secuencial aumentado MCMC 1: Iniciar con M valores iniciales (θm0 , x m 0 ) de la densidad a priori p(θ, x0). 2: Para t = 1, . . . , T hacer 3: Ajustar modelo p(θ | y0:t) con {1θBt−1, . . . , Mθ B t−1} 4: Para m = 1, . . . ,M hacer 5: Para k = 1, . . . , B hacer 6: Generar θ∗ ∼ q(θ∗ |θk−1) 7: Obtener x∗t ∼ p(x∗t |xt−1, θ ∗) 8: Calcular α = mı́n ( 1, p(yt |x∗t ,θ∗) p̂(θ∗ | y0:t−1) p(yt |xkt ,θk) p̂(θk | y0:t−1) ) 9: Hacer θk = { θ∗ con probabilidad mı́n(α, 1) θk−1 con probabilidad 1−mı́n(α, 1). 10: Fin Para 11: Fin Para 12: Recolectar {1θBt , . . . ,Mθ B t } 13: Generar xmt ∼ p(xt |xmt−1, θ m t ) 14: Fin Para Para los valores iniciales de x0 y y0 se usaron las distribuciones (3.9) y (3.10). Para la obtención de los pesos wt = p(yt |xt,θ) se utilizó la medida de probabilidad normal con promedio yt y desviación estándar 1, esto es N(· | yt, 1). Se usó esta probabilidad porque se obtuvieron mejores resultados que al usar la probabilidad dada por (3.12). Para obtener los valores iniciales de θ0 se usó la distribución normal multivaria- da truncada con valor mı́nimo (0, 0, 0, 0) y valor máximo (1, 1, 1, 1), con esperanza (0, 0, 0.75, 0), como matriz de covarianza una matriz diagonal con elementos 3.4. Algoritmos 35 (0.0005, 0.001, 0.05, 0.0000015). La distribución propuesta q que se utilizó fue una dis- tribución normal multivariada truncada con valor mı́nimo (0, 0, 0, 0) y valor máximo (1, 1, 1, 1), cuya matriz de covarianza es una matriz diagonal con elementos (0.06, 0.05, 0.05, 0.05). Se decidió usar estos valores pues al usar otros valores los resul- tados del ajuste no vaŕıan mucho. Este algoritmo permite realizar su implementación en paralelo. El ciclo que se realiza sobre el número de part́ıculas, cada iteración corresponde a la construcción de la cadena de Markov de la correspondiente part́ıcula m y es independiente de las otras. Se usó la función mclapply del paquete parallel1 de R. El tiempo de ejecución se redujo en un 80 % con respecto al tiempo de una implementación estándar de este mismo algoritmo. 3.4.5. Filtro de part́ıculas secuencial aumentado MCMC con remuestreo En esta versión del Algoritmo 3.3 filtro de part́ıculas secuencial aumentado MCMC se le agregó el paso de remuestreo de las muestras, con el propósito de mejorar la estimación de los caminos de la variable latente y la variable observada. Las ĺıneas del Algoritmo 3.4 se agrega al Algoritmo 3.3 después de la ĺınea 13. Algoritmo 3.4 Remuestreo en SAMCMC 1: Obtener los pesos ω (i) t ∼ p(yt |x(i) t ,θ (i) t ) 2: Realizar remuestreo con reemplazo de tamaño N del conjunto {x(i) t ,θ (i) t : i = 1, . . . , N} con probabilidades {ω(i) t : i = 1, . . . , N}. Para la obtención de los pesos se usa la medida de probabilidad (3.12). El remuestreo se realiza con el método usado en [19] llamado remuestreo sistemático, el cual tiene mejor varianza que el simple remuestreo aleatorio, además se espera que tenga un mejor 1Versión 3.6.3 del paquete parallel. 3.4. Algoritmos 36 desempeño que el remuestreo estratificado [16]. Es un remuestreo popular por su fácil implementación y por su baja complejidad computacional. El Algoritmo 3.5 muestra el procedimiento. Algoritmo 3.5 Remuestreo sistemático 1: Obtener u ∼ U([0, 1]) 2: Normalizar los pesos ω(i) 3: Hacer u = u/N , j = 1 y Sw = ω1. 4: Para K = 1, . . . , n hacer 5: Mientras Sw < u hacer 6: Hacer j = j + 1 y Sw = Sw + ωj. 7: Fin Mientras 8: Hacer ak = j y u = u+ 1/N . 9: Fin Para 10: Devolver a1:N Al remuestreo se le agrega una restricción debido a que el número de elementos únicos en la muestra de los parámetros se va reduciendo. La restricción consiste en aplicar el remuestreo si el número de elementos únicos en la muestra de los parámetros al realizar el remuestreo es mayor al 50 % del número total de part́ıculas. 37 CAPÍTULO 4 RESULTADOS Una vez que se contó con la información simulada y real, se procedió a aplicar los algoritmos y mostrar la eficacia en la estimación de los parámetros del modelo. El caṕıtulo se desarrolla de la siguiente manera: en la sección 4.1 se muestran los resultados del algoritmo filtro de part́ıculas Liu-West (LW), en la sección 4.2 los resultados del algoritmo filtro de part́ıculas Marginal Metropolis-Hastings (PMMH), la sección 4.3 contiene resultados del filtro Secuencial Aumentado MCMC (SAMCMC) además de pruebas de sensibilidad en el periodo de quema y número de part́ıculas, sección 4.4 se presentan los resultados del algoritmo SAMCMC agregándole remuestreo, en la sección 4.5 se comparan los resultados y se determina cual algoritmo es mejor y finalmente la sección 4.6 muestra los resultados de ejecutar SAMCMC en los datos reales. 4.1. Filtro de part́ıculas Liu-West El algoritmo de Liu-West (LW) es el punto de referencia para los otros dos algo- ritmos que se usaron en el estudio. El algoritmo usa 130 part́ıculas, sus resultados 4.1. Filtro de part́ıculas Liu-West 38 se presentan a continuación. En el Cuadro 4.1 se puede ver el resumen de valores estad́ısticos de la muestra obtenida. Cuadro 4.1: Resumen de valores estad́ısticos de las muestras obtenidas con LW. α β H µ Min. 0.1883 0.1911 0.4296 0.1760 5 % Per. 0.1889 0.1926 0.4303 0.1775 1st Qu. 0.1902 0.1965 0.4316 0.1815 Median 0.1909 0.2049 0.4323 0.1899 Mean 0.1922 0.2022 0.4337 0.1872 3rd Qu. 0.1934 0.2060 0.4351 0.1910 95 % Per. 0.1984 0.2080 0.4401 0.1930 Max. 0.1997 0.2085 0.4414 0.1935 El algoritmo generó muestras con poca variabilidad, pero la estimación de los parámetros del modelo no es aceptable, ni siquiera los valores reales se encuentran en el rango de las muestras, tal como esto se puede apreciar en los Cuadros 4.1 y 4.2. Cuadro 4.2: Valor promedio, desviación estándar y coeficiente de variación de la estimación de los parámetros junto con el valor real a estimar con LW. α β H µ Mean 0.19224 0.20215 0.43374 0.18715 Sd 0.00301 0.00531 0.00314 0.00533 CV 0.01566 0.02627 0.00724 0.02848 Real 0.02733 0.07567 0.60000 0.00140 4.1. Filtro de part́ıculas Liu-West 39 6.6 6.8 7.0 0 100 200 Día Lo g P re ci o Log Precio Stock, Camino simulado y estimación Figura 4.1: Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo LW A diferencia de la estimación de los parámetros, la estimación de la variable ob- servada Log Precio es muy buena. El algoritmo LW logró reconstruir muy bien el camino observado como se aprecia en la Figura 4.1 y se aprecia además que la banda de confianza a un 90 %, área sombreada color gris, es angosta. Se visualizan un par de observaciones que no entran en la banda pero aún aśı la estimación es aceptable. 4.1. Filtro de part́ıculas Liu-West 40 0.00 0.25 0.50 0.75 1.00 0 100 200 Día V ol at ili da d Volatilidad, Camino simulado y estimación Figura 4.2: Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja va- lor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Ĺıneas grises muestra de caminos estimados. Al- goritmo LW En cuanto a la estimación del camino de la volatilidad simulada, la cual se muestra en la Figura 4.2, vemos que el algoritmo no reconstruye muy bien la volatilidad obser- vada en los primeros 60 d́ıas aproximadamente pero después la estimación es buena. Tiene una banda de confianza que al inicio es muy ancha pero luego disminuye y man- tiene un ancho relativamente fijo. Por otro lado, el camino real se sale de la banda en pequeños tramos. Las medidas de la estimación, el RECM, IS90 % y IS80 % se presentan en el Cuadro 4.1. Filtro de part́ıculas Liu-West 41 4.3, más adelante se usarán para comparar con los otros algoritmos. Cuadro 4.3: Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 80 % y 90 % de las variables del modelo. Algoritmo LW. Variable RECM IS90 % IS80 % Volatilidad 0.16505 0.14641 0.11160 Log Precio 0.04263 1.42452 0.86378 α 0.10797 2.09150 1.05615 β 0.06824 1.27547 0.64866 H 0.27886 5.50899 2.76288 µ 0.13470 2.60580 1.31392 El tiempo promedio de las 10 ejecuciones es de 6.60 minutos. 4.1.1. Sensibilidad al número de part́ıculas El algoritmo no genera estimaciones aceptables de los parámetros del modelo, pero se realizaron escenarios variando el número de part́ıculas para ver que efectos tienen en los resultados. Los resultados se pueden observar en los Cuadros 4.4, 4.5. Cuadro 4.4: Coeficiente de variación de los escenarios de sensibilidad al número de part́ıcu- las en LW. Part́ıculas α β H µ 130 1.57 % 2.63 % 0.72 % 2.85 % 500 8.92 % 5.52 % 15.82 % 6.97 % 1,000 5.94 % 3.21 % 18.86 % 4.70 % 10,000 9.39 % 13.13 % 30.39 % 8.22 % Podemos ver que: 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 42 • En general es un algoritmo que consume poco tiempo en su ejecución. • El aumento en el número de part́ıculas no mejora las estimaciones. • Las muestras se vuelven más heterogéneas conforme aumenta el número de part́ıcu- las, es decir hay más variabilidad en las muestras. Cuadro 4.5: Valores promedio y desviación estándar para diferentes escenarios de número de part́ıculas en LW. Part́ıculas Tiempo (Minutos) α β H µ 130 6.75 mean 0.19224 0.20215 0.43374 0.18715 sd 0.00301 0.00531 0.00314 0.00533 500 25.69 mean 0.14205 0.14414 0.35632 0.13709 sd 0.01267 0.00795 0.05636 0.00955 1,000 52.32 mean 0.16584 0.17957 0.40489 0.16506 sd 0.00985 0.00577 0.07635 0.00776 10,000 528.83 mean 0.15171 0.14915 0.43330 0.14781 sd 0.01425 0.01958 0.13168 0.01215 Real 0.02733 0.07567 0.60000 0.00140 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings El algoritmo filtro de part́ıculas Marginal Metropolis-Hastings (PMMH) se ejecutó con un tamaño de muestra de 10,000 y se quemó1 el 50 %. Se usaron distribuciones normales para las distribuciones a priori y la propuesta. Como se puede apreciar en el Cuadro 4.6 y Figura 4.3, las muestras tienen histo- gramas asimétricos. Se ve la gran influencia de la distribución a priori, la cual tiene 1Periodo en el cual la cadena de Markov se mueve de su posición inicial a la región del espacio de parámetros que tiene alta probabilidad posterior. 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 43 una media (α, β,H, µ) = (0, 0, 0.75, 0). Se ve en la Figura 4.3 que los valores reales de (α, β, µ) se encuentran en zonas de alta probabilidad en las distribuciones emṕıri- cas de las muestras. Para el parámetro H su valor real se encuentra en una zona con probabilidad no tan alta en la distribución emṕırica de la muestra. Cuadro 4.6: Resumen de valores estad́ısticos de las muestras obtenidas con PMMH. α β H µ Min. 0.000001 0.0000654 0.3845 0.0000007 5 % Per. 0.00446 0.0111764 0.5859 0.000128 1st Qu. 0.01465 0.0439644 0.6869 0.000520 Median 0.02746 0.0765881 0.7496 0.000979 Mean 0.03116 0.0884717 0.7461 0.001143 3rd Qu. 0.04513 0.1177788 0.8102 0.001525 95 % Per. 0.06824 0.2074704 0.8952 0.002806 Max. 0.08825 0.3571898 0.9951 0.004197 El valor estimado de los parámetros, desviación estándar y coeficiente de variación se muestra en el Cuadro 4.7. Se obtienen estimaciones buenas, a excepción del parámetro H. Con respecto a la variabilidad, en general es alta. Cuadro 4.7: Valor promedio, desviación estándar y coeficiente de variación de la estimación de los parámetros junto con el valor real a estimar con PMMH. α β H µ Mean 0.03116 0.08847 0.74609 0.00114 Sd 0.01998 0.06209 0.09437 0.00082 CV 0.64121 0.70182 0.12649 0.71930 Real 0.02733 0.07567 0.60000 0.00140 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 44 0 5 10 15 20 0.000 0.025 0.050 0.075 α D en si da d Densidad de Alfa 0 2 4 6 8 0.0 0.1 0.2 0.3 β D en si da d Densidad de Beta 0 2 4 0.4 0.6 0.8 1.0 H D en si da d Densidad de H 0 200 400 600 0.000 0.001 0.002 0.003 0.004 µ D en si da d Densidad de Mu Figura 4.3: Histogramas de las muestras obtenidas con PMMH. La ĺınea discontinua de color rojo muestra el valor real. El área sombreada color gris corresponde a la estimación de la densidad por medio de un kernel gaussiano. 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 45 La Figura 4.3 une la información de los dos cuadros anteriores. Exhibe los histo- gramas y la área sombreada que permite visualizar la densidad emṕırica de las dis- tribuciones. Una ĺınea roja vertical muestra donde se ubica el valor real dentro de la densidad. No todos los valores reales se ubican en zonas de alta probabilidad posterior, algo esperado en una estimación aceptable. 0.000 0.025 0.050 0.075 5000 6000 7000 8000 9000 10000 Iteración α Traza de Alfa 0.0 0.1 0.2 0.3 5000 6000 7000 8000 9000 10000 Iteración β Traza de Beta 0.4 0.6 0.8 1.0 5000 6000 7000 8000 9000 10000 Iteración H Traza de H 0.000 0.001 0.002 0.003 0.004 5000 6000 7000 8000 9000 10000 Iteración µ Traza de Mu Figura 4.4: Gráfico de traza de las muestras obtenidas con PMMH. La ĺınea discontinua de color rojo muestra el valor real. Los gráficos de traza del MCMC2 para cada parámetro se ve en la Figura 4.4. La traza de α da evidencia de tener ciertas tendencias y aunque se quemaron las 2Siglas del inglés Markov Chain Monte Carlo. 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 46 primeras 5,000 iteraciones, no parece tener una buena mezcla. Los procesos de los otros parámetros tienen en general buena mezcla. 6.6 6.7 6.8 6.9 7.0 7.1 0 100 200 Día Lo g P re ci o Log Precio Stock, Camino simulado y estimación Figura 4.5: Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo PMMH. En la Figura 4.5 se aprecia que la reconstrucción del camino del Log Precio no es buena, pues hay varias secciones del camino que quedan por fuera de la banda de 90 % de confianza. De forma adicional, se observa que prácticamente no hay diferencia entre usar el promedio o la mediana como estimador. 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 47 0.0 0.2 0.4 0.6 0 100 200 Día V ol at ili da d Volatilidad, Camino simulado y estimación Figura 4.6: Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja valor promedio de la estimación, banda gris intervalo de confianza 90 %. Algoritmo PMMH. Al igual que en los otros algoritmos, la estimación del camino de la volatilidad no es buena. Como se aprecia en la Figura 4.6 la estimación decrece con el tiempo, situación contraria al camino real. La banda de confianza es ancha prácticamente desde el inicio hasta el final del camino. 4.2. Filtro de part́ıculas Marginal Metropolis-Hastings 48 Cuadro 4.8: Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 80 % y 90 % de las variables del modelo. Algoritmo PMMH. Variable RECM IS90 % IS80 % Volatilidad 0.20129 0.55102 0.41654 Log Precio 0.10597 0.47259 0.37931 α 0.01777 0.05471 0.04313 β 0.06295 0.19630 0.15490 H 0.17194 0.38390 0.52944 µ 0.00082 0.00234 0.00183 Las medidas del desempeño del algoritmo se muestran en el Cuadro 4.8. Se vi- sualizan valores altos para la volatilidad, log Precio y H. En la sección 4.5 se verá el desempeño en relación a los otros algoritmos. El tiempo promedio de las 10 ejecuciones es de 11.84 horas. 4.2.1. Distribución a priori Gamma Se realizó un escenario en el cual se cambió la distribución conjunta a priori de los parámetros, con el fin de ver si la influencia en la estimación se sigue presentando. Se es- cogió una distribución gamma con con parámetro de forma (α, β,H, µ) = (5, 13, 100, 2) y parámetro de tasa (α, β,H, µ) = (150, 150, 150, 1500). Con estos parámetros la es- peranza y la desviación estándar son (α, β,H, µ) = (0.0333, 0.0867, 0.6667, 0.0013) y (α, β,H, µ) = (0.0149, 0.0240, 0.0667, 0.0009) respectivamente. Como se puede ver en los Cuadros 4.9 y 4.10, el promedio y la desviación estándar de las muestras de los parámetros prácticamente son los de la distribución a priori. Esto se puede ver como una debilidad del algoritmo, pues si no se cuenta con información a priori se tiene que hacer uso de distribuciones no informativas [30] que no van a ayudar a obtener buenas estimaciones. 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 49 Cuadro 4.9: Resumen de valores estad́ısticos de las muestras obtenidas con PMMH Gamma. α β H µ Min. 0.005424 0.03028 0.4700 0.000011 5 % Per. 0.016273 0.05348 0.5551 0.000239 1st Qu. 0.023986 0.07031 0.6184 0.000688 Median 0.034849 0.08375 0.6589 0.001268 Mean 0.036093 0.08666 0.6571 0.001566 3rd Qu. 0.046388 0.10153 0.6974 0.002351 95 % Per. 0.061761 0.12619 0.7559 0.003699 Max. 0.073953 0.16053 0.8859 0.005089 Cuadro 4.10: Valor promedio, desviación estándar y coeficiente de variación de la estima- ción de los parámetros junto con el valor real a estimar con PMMH Gamma. α β H µ Mean 0.03609 0.08666 0.65711 0.00156 Sd 0.01454 0.02261 0.06007 0.00109 CV 0.40288 0.26090 0.09142 0.69872 Real 0.02733 0.07567 0.60000 0.00140 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC En esta sección se presenta los resultados de aplicar el algoritmo filtro de part́ıculas Secuencial Aumentado MCMC (SAMCMC) en los datos simulados. El algoritmo usa 130 part́ıculas que es equivalente al número de MCMC’s que se ejecutan en cada iteración de t. Para cada MCMC que realiza un periodo de quema de 1,000 iteraciones. 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 50 0 30 60 90 0.03 0.04 0.05 α D en si da d Densidad de Alfa 0 10 20 0.05 0.10 0.15 β D en si da d Densidad de Beta 0 3 6 9 12 0.50 0.55 0.60 0.65 0.70 H D en si da d Densidad de H 0 500 1000 1500 0.001 0.002 µ D en si da d Densidad de Mu Figura 4.7: Histogramas de las muestras obtenidas con SAMCMC. La ĺınea discontinua de color rojo muestra el valor real. El área sombreada color gris corresponde a la estimación de la densidad por medio de un kernel gaussiano. 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 51 El principal resultado del algoritmo es una muestra de la distribución posterior de cada uno de los parámetros del modelo de espacio estado (2.11), (2.12). En el Cuadro 4.11 se muestra un resumen de valores estad́ısticos de la muestra de cada uno de los 4 parámetros. El cuadro da evidencia de que las distribuciones son simétricas. Cuadro 4.11: Resumen de valores estad́ısticos de las muestras obtenidas con SAMCMC. α β H µ Min.: 0.02216 0.01682 0.4972 0.0004258 5 % Pe.: 0.02470 0.05347 0.5349 0.0006655 1st Qu.: 0.02898 0.08060 0.5729 0.0010783 Median: 0.03259 0.09515 0.6109 0.0013529 Mean: 0.03299 0.09884 0.6103 0.0013438 3rd Qu.: 0.03593 0.11969 0.6457 0.0016134 95 % Pe.: 0.04326 0.14456 0.6868 0.0019719 Max.: 0.04936 0.18402 0.7253 0.0026080 Cuadro 4.12: Valor promedio y desviación estándar de la estimación de los parámetros junto con el valor real a estimar con SAMCMC. α β H µ Mean 0.03299 0.09884 0.61030 0.00134 Sd 0.00562 0.02888 0.04750 0.00042 CV 0.17035 0.29219 0.07783 0.31343 Real 0.02733 0.07567 0.6000 0.00140 Como los datos son simulados, se conocen los valores reales de los parámetros que los generaron. El Cuadro 4.12 presenta el valor estimado (el valor promedio de la muestra), la desviación estándar y el coeficiente de variación para los parámetros. Los valores 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 52 estimados son aceptables, el parámetro β es el que está un poco lejos del valor real y el parámetro µ es el que tiene mayor variabilidad con respecto a los otros parámetros. La Figura 4.7 se aprecia que los valores reales se ubican en zonas de alta probabilidad posterior, algo esperado en una estimación aceptable. 5.5 6.0 6.5 7.0 7.5 0 100 200 Día Lo g P re ci o Log Precio Stock, Camino simulado y estimación Figura 4.8: Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo SAMCMC. El objetivo principal es la estimación de los parámetros del modelo, pero también se puede ver como el algoritmo se desempeñó reconstruyendo el camino del logaritmo de los precios. En la Figura 4.8, la ĺınea negra es el camino simulado del logaritmo de los 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 53 precios. La ĺınea rojo discontinua es el camino de los promedios para cada d́ıa, esto pues para cada d́ıa hay una muestra de la distribución posterior de la variable observada. La banda gris representa un intervalo de confianza de 90 %, que se obtiene con los intervalos diarios. La banda de confianza tiende a ensancharse conforme aumenta el tiempo. 0.00 0.25 0.50 0.75 1.00 0 100 200 Día V ol at ili da d Volatilidad, Camino simulado y estimación Figura 4.9: Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja valor promedio de la estimación, banda gris intervalo de confianza 90 %. Ĺıneas grises son muestras de caminos estimados. Algoritmo SAMCMC. Para el algoritmo lo más dif́ıcil es reconstruir el camino de la variable latente, es decir el camino de la volatilidad. La Figura 4.9 replica la idea de la Figura anterior 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 54 con la adición de 20 caminos generados por el algoritmo, se aprecia que hay mucha dispersión en los caminos. La banda de confianza es grande y casi con el mismo ancho en todo el periodo. A diferencia de los algoritmos LW y PMHM, este algoritmo no incorpora el paso de selección. En este paso se calcula la probabilidad p(yobt |x̂it, θ̂ i ) para cada part́ıcula y se realiza un remuestreo con reemplazo, con lo cual las part́ıculas más representativas quedan en la muestra del d́ıa t y se usan en la iteración siguiente t+ 1. Esta selección afecta directamente en la reconstrucción de los caminos, como se evidencia en las Figuras. El algoritmo LW, como se puede ver en el Algoritmo 3.1, realiza un remuestreo adicional. Este remuestreo se realiza al inicio de la iteración con el propósito de escoger las part́ıculas con mayor probabilidad de que cuando evolucionen en el tiempo t+1 sean más consistente con el valor observado yt+1 que las otras part́ıculas. Esta caracteŕıstica junto con el factor del error o aleatoriedad produce que las estimaciones de los caminos sean mejores a pesar de la deficiente estimación de los parámetros. Cuadro 4.13: Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 90 % y 80 % de las variables del modelo. Algoritmo SAMCMC. Variable RECM IS90 % IS80 % Volatilidad 0.31220 1.29910 0.92296 Log Precio 0.42067 0.89250 0.78645 α 0.01455 0.04371 0.03802 β 0.03224 0.09855 0.08763 H 0.08886 0.25959 0.20863 µ 0.00050 0.00140 0.00130 También hay que considerar que sobre esta variable latente recae toda la dispersión del sistema. Hay tres fuentes de dispersión (α, β,H), que son los parámetros de la 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 55 ecuación de estado. Las medidas a usar para las comparaciones entre algoritmos se muestran en el Cuadro 4.13. El tiempo promedio de las 10 ejecuciones es de 35.80 horas. 4.3.1. Sensibilidad al periodo de quema en los MCMC’s Se ejecutó el algoritmo varias veces cambiando el periodo de quema con el fin de ver su comportamiento. Los periodos de quema que se probaron fueron: 100, 500, 1500, 1000 (escenario base). El número de part́ıculas que se usó en cada escenario es de 130. Cuadro 4.14: Coeficiente de variación de los escenarios de sensibilidad al periodo de quema algoritmo SAMCMC. Quema α β H µ 100 43.07 % 43.87 % 19.84 % 49.45 % 500 47.22 % 39.73 % 14.16 % 35.34 % 1,000 17.04 % 29.22 % 7.78 % 31.34 % 1,500 25.42 % 31.12 % 11.22 % 25.40 % 2,000 29.09 % 21.25 % 13.92 % 36.90 % Los Cuadros 4.14 y 4.15 muestran los resultados de los parámetros aśı como el tiempo de ejecución. Se nota que: • El crecimiento del tiempo es lineal con respecto al periodo de quema. • Los parámetros α y β no tienen tan clara su tendencia, pero en general su des- viación estándar tiende a disminuir. • El valor promedio del parámetro H tiende a disminuir con el aumento del periodo de quema y la menor desviación estándar se obtuvo con 1,000. • El valor promedio del parámetro µ tiende a aumentar conforme aumenta el pe- riodo de quema y su desviación estándar es muy simular entre los escenarios. 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 56 • La muestra más homogénea se logró con un periodo de quema de 1,000. Cuadro 4.15: Valores promedio y desviación estándar para diferentes escenarios de periodo de quema en los MCMC’s algoritmo SAMCMC. Quema Tiempo (Horas) α β H µ 100 3.73 mean 0.02626 0.04169 0.64018 0.00091 sd 0.01131 0.01829 0.12701 0.00045 500 18.09 mean 0.02810 0.05178 0.62088 0.00116 sd 0.01327 0.02057 0.08794 0.00041 1,000 35.66 mean 0.03299 0.09884 0.61025 0.00134 sd 0.00562 0.02888 0.04750 0.00042 1,500 53.28 mean 0.03128 0.05354 0.55786 0.00189 sd 0.00795 0.01666 0.06259 0.00048 2,000 71.33 mean 0.02025 0.05530 0.49280 0.00168 sd 0.00589 0.01175 0.06862 0.00062 Real 0.02733 0.07567 0.60000 0.00140 4.3.2. Sensibilidad al periodo de quema y número de part́ıcu- las En estos escenarios se varió el número de part́ıculas entre 130, 500 y 1000, con un periodo de quema de 100. Aśı como número de part́ıculas 130, 150 y 1000 con periodo de quema de 1000. 4.3. Filtro de part́ıculas Secuencial Aumentado MCMC 57 Cuadro 4.16: Coeficiente de variación de los escenarios de sensibilidad al periodo de quema y número de part́ıculas. Algoritmo SAMCMC. Quema Part́ıculas α β H µ 100 130 43 % 44 % 20 % 49 % 100 500 48 % 49 % 23 % 45 % 100 1,000 44 % 43 % 21 % 50 % 1,000 130 17 % 29 % 8 % 31 % 1,000 150 23 % 24 % 12 % 29 % 1,000 1,000 32 % 28 % 17 % 28 % Cuadro 4.17: Valores promedio y desviación estándar para diferentes escenarios de periodo de quema y número de part́ıculas en los MCMC’s. Algoritmo SAMCMC. Quema Part́ıculas Tiempo (Horas) α β H µ 100 130 3.73 mean 0.02626 0.04169 0.64018 0.00091 sd 0.01131 0.01829 0.12701 0.00045 100 500 13.71 mean 0.02943 0.03461 0.63181 0.00129 sd 0.01411 0.01683 0.14370 0.00058 100 1,000 27.33 mean 0.02533 0.03679 0.65284 0.00133 sd 0.01122 0.01579 0.13501 0.00066 1,000 130 35.66 mean 0.03299 0.09884 0.61025 0.00134 sd 0.00562 0.02888 0.04750 0.00042 1,000 150 40.13 mean 0.04528 0.08410 0.61254 0.00113 sd 0.01048 0.02001 0.07613 0.00033 1,000 1,000 263.28 mean 0.03225 0.05652 0.57516 0.00139 sd 0.01029 0.01593 0.09794 0.00039 Real 0.02733 0.07567 0.60000 0.00140 En los Cuadros 4.17 y 4.16 se muestras los resultados de los escenarios. Podemos 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo 58 observar que: • El crecimiento del tiempo es lineal con respecto al número de part́ıculas. • Como era de esperar, el aumento en el periodo de quema produce muestras más homogéneas. • Para los parámetros α y H con un periodo de quema de 1000, se observa que al aumentar el número de part́ıculas las muestras tienden a aumentar la variabilidad. • En el caso de β y µ con periodo de quema de 1000, el aumento en el número de part́ıculas no muestra un tendencia clara. 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo Dado que el algoritmo SAMCMC se implementó sin el paso de remuestreo y se obtuvo estimaciones buenas para los parámetros pero malas para la estimación de los caminos se procedió a implementar el algoritmo con el paso de remuestreo. Un resumen de los valores estad́ısticos se muestran en el Cuadro 4.18. En la Figura 4.10 se presenta los histogramas de las muestras obtenidas. El remuestreo provoca una reducción de elementos únicos en la muestra esto se puede apreciar al comparar la Figura 4.7 y 4.10. En comparación con SAMCMC el agregar remuestreo desmejora la estimación de los parámetros, ver Cuadro 4.19. 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo 59 0 50 100 150 200 0.01 0.02 α D en si da d Densidad de Alfa 0 20 40 60 0.000 0.025 0.050 0.075 0.100 0.125 β D en si da d Densidad de Beta 0 3 6 9 0.5 0.6 0.7 0.8 0.9 H D en si da d Densidad de H 0 1000 2000 3000 4000 0.0005 0.0010 0.0015 µ D en si da d Densidad de Mu Figura 4.10: Histogramas de las muestras obtenidas con SAMCMC con remuestreo. La ĺınea discontinua de color rojo muestra el valor real. El área sombreada color gris corresponde a la estimación de la densidad por medio de un kernel gaussiano. 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo 60 Cuadro 4.18: Resumen de valores estad́ısticos de las muestras obtenidas con SAMCMC con remuestreo. α β H µ Min.: 0.004734 0.003275 0.4401 0.0002320 5 % Pe.: 0.007247 0.025687 0.4971 0.0007125 1st Qu.: 0.011274 0.048534 0.6171 0.0009292 Median: 0.014037 0.068791 0.6733 0.0010446 Mean: 0.014167 0.064153 0.6735 0.0011197 3rd Qu.: 0.017026 0.081621 0.7381 0.0013254 95 % Pe.: 0.020627 0.094687 0.7987 0.0016118 Max.: 0.023010 0.120786 0.8946 0.0018735 Cuadro 4.19: Valor promedio y desviación estándar de la estimación de los parámetros junto con el valor real a estimar con SAMCMC con remuestreo. α β H µ Mean 0.01417 0.06415 0.6735 0.00112 Sd 0.00412 0.02223 0.0949 0.00029 CV 0.29076 0.34653 0.1409 0.25893 Real 0.02733 0.07567 0.6000 0.00140 El beneficio que se obtiene por el remuestreo es en la estimación del camino del logaritmo del precio y el camino de la estimación de la volatilidad. Esto se puede ver en las Figuras 4.11 y 4.12. 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo 61 6.6 6.8 7.0 0 100 200 Día Lo g P re ci o Log Precio Stock, Camino simulado y estimación Figura 4.11: Camino del Log Precio Stock simulado. Ĺınea discontinua roja valor promedio de la estimación. Ĺınea discontinua verde es mediana, banda gris intervalo de confianza 90 %. Algoritmo SAMCMC con remuestreo. 4.4. Filtro de part́ıculas Secuencial Aumentado MCMC con remuestreo 62 0.25 0.50 0.75 0 100 200 Día V ol at ili da d Volatilidad, Camino simulado y estimación Figura 4.12: Camino de la volatilidad simulada ĺınea color azul. Ĺınea discontinua roja va- lor promedio de la estimación, banda gris intervalo de confianza 90 %. Ĺıneas grises son muestras de caminos estimados. Algoritmo SAMCMC con remues- treo. Las métricas respectivas para este algoritmo se muestran en el Cuadro 4.20. El tiempo promedio de las 10 ejecuciones es de 35.72 horas. 4.5. Resumen de resultados 63 Cuadro 4.20: Ráız del error cuadrático medio (RECM) e Interval score (IS) con niveles 90 % y 80 % de las variables del modelo. Algoritmo SAMCMC con remuestreo. Variable RECM IS90 % IS80 % Volatilidad 0.14114 0.39840 0.30789 Log Precio 0.09486 0.65224 0.51636 α 0.01287 0.10432 0.06461 β 0.03756 0.43383 0.25860 H 0.11488 0.28125 0.26501 µ 0.00059 0.00162 0.00140 4.5. Resumen de resultados En esta sección se comparan el tiempo de ejecución, variabilidad, la ráız del error cuadrático medio y el score de intervalo para los diferentes algoritmos. El algoritmo que dura menos es LW con 6.60 minutos, esto se debe a que el algoritmo sólo itera sobre el tiempo t = 1, . . . , 255. El algoritmo PMMH tiene un tiempo de ejecución de 11.84 horas, cabe recalcar que este algoritmo ejecuta el filtro de part́ıculas estándar 10,000 veces. El método SAMCMC tiene un tiempo de ejecución de 35.80 horas y SAMCMC con remuestreo un tiempo de 35.72 horas, son los más exigentes en cálculos. Estos últimos realizan en cada iteración t de la variable tiempo 130 MCMC, es decir una cadena para cada part́ıcula, cada una con 1,000 iteraciones de periodo de quema. En cuanto a la variabilidad de la estimación de los parámetros, el algoritmo de LW tiene los coeficientes de variación más bajas, pero la estimación que se obtiene es la más deficiente de todas. Con respecto a los otros algoritmos es SAMCMC el que tiene menor variabilidad. Por ejemplo el parámetro H en el algoritmo SAMCMC se obtuvo un coeficiente de variación de 0.07783, mientras que en los algoritmos PMMH y 4.5. Resumen de resultados 64 SAMCMC con remuestreo prácticamente se duplico la variabilidad, pasando a 0.12649 y 0.1409 respectivamente. Cuadro 4.21: Comparación de la ráız del error cuadrático medio. Promedios sobre las 10 ejecuciones. RECM LW PMMH SAMCMC SAMCMCR Volatilidad 0.1650 0.2013 0.3122 0.1411 Log Precio 0.0426 0.1060 0.4207 0.0949 α 0.1080 0.0178 0.0146 0.0129 β 0.0682 0.0630 0.0322 0.0376 H 0.2789 0.1719 0.0889 0.1149 µ 0.1347 0.0008 0.0005 0.0006 Con respecto a la ráız del error cuadrático medio, ver Cuadro 4.21, vemos que la volatilidad se estimó mejor con SAMCMC con remuestreo y el Log precio se estimó mejor con LW pero como se indicó su estimación de parámetros es deficiente. La se- gunda mejor estimación es LW para la volatilidad y SAMCMC con remuestreo para el Log Precio. Por otra parte, es SAMCMC en general el que tiene mejor RECM en la estimación de los parámetros y SAMCMC con remuestreo el segundo mejor. En los Cuadros 4.22 y 4.23 se ve que SAMCMC es el algoritmo que tiene mejores valores de IS90 % y IS80 % para los parámetros. Estos valores indican que su banda de confianza es mejor predictora en relación al resto de algoritmos. Y en general es SAMCMC con remuestreo el que tiene el segundo lugar en mejores valores en todos los variables. En resumen, conforme se aumenta la complejidad de los algoritmos las estimaciones de los caminos de las variables Log precio y la volatilidad desmejoran, mientras que la estimación de los parámetros del modelo mejora. Es el algoritmo SAMCMC con 4.5. Resumen de resultados 65 remuestreo el que tiene un balance de medidas buenas en todas las variables. Cuadro 4.22: Comparación del score de intervalo. Promedios sobre las 10 ejecuciones. IS90 % LW PMMH SAMCMC SAMCMCR Volatilidad 0.1464 0.5510 1.2991 0.3984 Log Precio 1.4245 0.4726 0.8925 0.6522 α 2.0915 0.0547 0.0437 0.1043 β 1.2755 0.1963 0.0986 0.4338 H 5.5090 0.3839 0.2596 0.2813 µ 2.6058 0.0023 0.0014 0.0016 Cuadro 4.23: Comparación del score de intervalo. Promedios sobre las 10 ejecuciones. IS80 % LW PMMH SAMCMC SAMCMCR Volatilidad 0.1116 0.4165 0.9230 0.3079 Log Precio 0.8638 0.3793 0.7865 0.5164 α 1.0562 0.0431 0.0380 0.0646 β 0.6487 0.1549 0.0876 0.2586 H 2.7629 0.5294 0.2086 0.2650 µ 1.3139 0.0018 0.0013 0.0014 En conclusión, el mejor algoritmo para estimar los parámetros del modelo es SAMCMC, pues tiene los valores más bajos en RECM, IS y variabilidad. Pero es SAMCMC con remuestreo en general el mejor algoritmo porque logra estimar mejor los caminos y parámetros. 4.6. Aplicación con datos reales 66 4.6. Aplicación con datos reales En esta parte del estudio se tomó el mejor algoritmo, Secuencial Aumentado MCMC (SAMCMC) con remuestreo, y se aplicó en una serie corta y otra larga del ı́ndice S&P500. Se usó este ı́ndice porque el art́ıculo [6] indica que hay evidencia de memoria larga, además porque usaron los datos de la serie corta para aplicar el procedimiento presentado en el art́ıculo. Los datos se encuentran en el anexo A.2. 4.6.1. Datos de un mes La serie corta del ı́ndice S&P500 corresponde a 22 observaciones diarias entre las fechas 2008-12-29 y 2009-01-29. El Cuadro 4.24 presenta los valores promedios estima- dos, 3 percentiles, la desviación estándar y el coeficiente de variación. El parámetro H presenta la variabilidad más baja al tener un coeficiente de variación de 0.20, mientras que los otros parámetros su variabilidad es mayor pero similar entre ellos. Cuadro 4.24: Valores estad́ısticos de los parámetros obtenidos con SAMCMC con remues- treo en datos de S&P500 del periodo 2008-12-29 al 2009-01-29. Los valores de Chronopoulou corresponden al estudio realizado en [6]. α β H µ mean 0.01896 0.04585 0.69483 0.00114 sd 0.00925 0.02118 0.13986 0.00047 CV 0.48787 0.46194 0.20129 0.41228 5 % Pc 0.00645 0.01448 0.43830 0.00039 50 % Pc 0.01966 0.04567 0.70971 0.00117 95 % Pc 0.03239 0.08142 0.94383 0.00186 Chronopoulou 0.05850 NA 0.53000 0.00540 Con esta misma serie de datos el art́ıculo de [6] estimó valores para α, H y µ, los cuales se muestran en el Cuadro 4.24. A diferencia del art́ıculo [6] en este trabajo se 4.6. Aplicación con datos reales 67 estiman los parámetros de forma conjunta. El tiempo que tomó la ejecución fue de 35.28 minutos. La Figura 4.13 presenta el camino del logaritmo del ı́ndice S&P500 junto a la esti- mación que corresponde al valor promedio de las muestras obtenidas por el algoritmo para cada uno de los d́ıas. La banda gris muestra una zona de confianza de 90 %. Se tiene una buena estimación. En la Figura 4.14 se presenta el camino de la volatilidad, construida con los pro- medios de las muestras diarias, además se presenta la mediana. 6.70 6.75 6.80 6.85 6.90 Dec 29 Jan 05 Jan 12 Jan 19 Jan 26 Fecha Lo g P re ci o Log Precio Stock Figura 4.13: Camino del Log S&P500 del periodo 2008-12-29 al 2009-01-29. Ĺınea disconti- nua roja valor promedio de la estimación. La ĺınea discontinua verde mediana, banda gris intervalo de confianza 90 %. 4.6. Aplicación con datos reales 68 0.00 0.25 0.50 0.75 0 5 10 15 20 Días V ol at ili da d Volatilidad estimada Figura 4.14: Camino de la volatilidad estimada usando los datos de S&P500 del periodo 2008-12-29 al 2009-01-29. Ĺınea discontinua roja valor promedio y la ĺınea discontinua verde mediana. 4.6.2. Datos de un año Para la serie larga del ı́ndice S&P500 se tomó información entre las fechas 2008- 12-29 y 2009-12-28, en total son 252 observaciones. El Cuadro 4.25 muestra los valores promedios estimados, 3 percentiles, la desviación estándar y el coeficiente de variación. Las desviaciones estándar disminuye: por ejemplo, el parámetro H tuvo una dismi- nución de 64 % con respecto a la estimació