Estimación puntual de XARIMA con prior y restricciones y mas ...

[Jorge]

Se solicita la implementación de método de estimación puntual, posiblemente basado en nlopt de modelos con las siguientes características:

• estructuras ARIMA multiestacionales
• inputs lineales: parámetros beta
• restricciones lineales sobre los parámetros beta
• prior conocido sobre los parámetros beta
• omitidos en el output (sin prior y devueltos en forma de serie como lo hace Estimate)
• matriz de covarianzas asintótica aproximada mediante el jacobiano s^2 (J'J)^-1

En una segunda fase:

• efectos no lineales delta
• efectos aditivos en terminos originales
• omitidos en los inputs
• prior de omitidos

[Víctor de Buen Remiro]

La solución a este problema se enmarcará dentro del paquete MPM (Modular Probabilistic Models) que permite definir de forma modular modelos estadísticos de cualquier tipo, ensamblando componentes para poder reutilizar el código al máximo.

Los diferentes módulos se irán añadiendo conforme se vayan necesitando, así como los métodos de estimación máximo-verosímiles y bayesianos más adecuados a cada caso.

Por el momento ya consta de estos estimadores genéricos

• método de optimización SBPLX de NLOPT del paquete NonLinGloOpt que no requiere el cálculo del gradiente y que admite restricciones lineales y no lineales.
• método Random Walk Metropolis Hastings del paquete BysSampler

Los módulos disponibles ahora mismo son

• bloque lineal con restricciones de rango
• bloque ARIMA
• bloque de varianza

He hecho una batería de pruebas de modelos X-ARIMA generados arbitrariamente por ingeniería inversa para comparar resultados de MPM+NLOPT con Estimate. Los resultados son siempre favorables a MPM+NLOPT en términos de verosimilitud. Los timepos de ejecución son más bajos en Estimate para modelos pequeños (menos de 100 variables), debido al overhead del manejo de estructuras de datos complejas en TOL, pero con medelos pequeños y grandes Estimate tarda bastante más, pues ahí la CPU está dominada por las operaciones matriciales y aritméticas.

También he probado un modelo real de distribución de prensa que es bastante grande con unas 500 variables y 4000 datos y MPM+NLOPT tarda la mitad (508 seg.) que Estimate (1042 seg.). Los resiudos son casi los mismos en ambos casos pero las soluciones no se parecen en nada, y eso pasa cuando hay mucha multicolinealidad.

Al introducir las restricciones de rango en ese problema real con MPM+NLOPT pasa de tardar unos minutos a 1 hora y 24 minutos. El error sube de 0.0527 a 0.0591, que puede ser razonable, pero hay muchas variables que se sitúan en la frontera.

No sé muy bien a qué puede deberse este problema, pero se me ocurren varias causas posibles:

• si es que está mal diseñado el modelo y las restricciones no son compatibles con los datos,
• si es que no se alcanza la convergencia por algún error de implementación mío o de NLOPT, o
• si simplemente los métodos libres de gradiente no son adecuados para este tipo de modelos con restricciones y multicolinealidad, que son dos cosas que hacen siempre bastante mala mezcla.

El gradiente de los modelos ARIMA y derivados, al igual que la mayoría de los modelos medianamente complejos, no es calculable analíticamente y en otros casos en que pudiera hacerse muy probablemente sería demasiado costoso. Debido a las restricciones de estacionariedad y otras que pudieran incorporarse, tampoco es trivial el cálculo del gradiente numérico aunque es algo que voy a intentar incorporar para poder hacer uso de métodos diferenciales de optimización como los que aportan NLOPT e IPOPT.

Otra posibilidad que merece la pena probar es el paquete ​levmar (Levenberg-Marquardt nonlinear least squares algorithms in C/C++) que permite usar métodos de diferencias finitas para el cálculo del Jacobiano. Este método sería válido para cualquier modelo de regresión, pero no para otros tipos de modelos más generales.

[Víctor de Buen Remiro]

Los problemas de convergencia se atenúan bastante mediante la introducción de información a priori complementaria para las restricciones de intervalo acotadas por ambos lados.

En este caso se trataba siempre de promociones para las que se introduce un límite superior arbitrariamente situado en 1 aunque en realidad es casi imposible que una promoción alcance siquiera un coeficiente de 80%. El límite inferior es cero porque en principio nadie va a dejar de comprar el periódico porque regalen algo. Si añadimos un prior normal con media 15% y desviación 20% entonces la densidad del 80% es insignificante:

El efecto de estos priors es evitar que la función de verosimilitud sea tan plana como es a causa de la multicolinealidad y la escasa superficie de respuesta que está entorno a 8, pues aunque hay unos 4000 datos hay unas 500 variables contando los omitidos.

Para que el modelo fuera más consistente habría que usar estructuras de tipo jerárquico para hacer depender el efecto de cada promoción de criterios más objetivos, como el valor o el interés del regalo y las condiciones para su obtención por parte del lector.

En general habría que dotar al usuario de mecanismos de diseño de información a priori prefabricada para las situaciones más habituales y enseñar a los analistas avanzados a crear sus propios diseños.

En cualquier caso la estimación puntual sigue topándose en algunos casos con las fronteras y sólo nos informa de la moda, que en modelos con tantas restricciones nos aporta muy poca información acerca de la distribución.

Es indispensable por tanto el uso de estimación bayesiana MCMC. Para ello se ha implementado un generador de candidatos asimétrico para el método Metropolis-Hastings que muestrea normales truncadas entre los límites de cada coordenada que ha sido previamente escalada para mejorar el mezclado de la cadena.