La simulación de parámetros y valores pronosticados a partir de la distribución a posteriori y predictiva se realiza con técnicas como Monte Carlo por Cadenas de Markov, introducida en la siguiente sección.

Monte Carlo por Cadenas de Markov (MCMC)

El método de Monte Carlo genera valores independientes tomando como base una distribución de probabilidad deseada. Las cadenas de Markov se basan en una estructura de dependencia entre los valores simulados consecutivamente. Al unir ambas técnicas se crea el método Monte Carlo por Cadenas de Markov, que consiste en hacer un muestreo a partir de distribuciones de probabilidad basadas en la construcción de cadenas, donde cada valor simulado, de acuerdo a [12] y [16], tiene dependencia con el dato anterior, llegando a una convergencia a la distribución deseada. Después de gran cantidad de simulaciones se establecen resultaos que se utilizan como una muestra incorrelacionada de la distribución deseada.

Los métodos MCMC son una estrategia para generar valores de la variable aleatoria buscada X mientras explora el espacio usando cadenas de Markov.

La simulación para actualizar el valor del parámetro X para este trabajo se llevó a cabo con el proceso de MCMC, que se basa en el muestreador de Gibbs, y este, a su vez, utiliza las distribuciones a posteriori condicionales conocidas de los parámetros, valor que permite la actualización de la función predictiva para el pronóstico.

Muestreador de Gibbs

Es una de las técnicas más usadas para las simulaciones de procesos de Monte Carlo por Cadenas de Markov (MCMC). El muestreador de Gibbs requiere conocimiento específico sobre la naturaleza condicional de la relación entre variables de interés. La idea básica es realizar un muestreo para el parámetro, por ejemplo, X, a partir de su distribución a posteriori condicional. Cuando se involucran más parámetros, este muestreador utiliza cada una de las distribuciones a posteriori condicionales por parámetro.

Este método se puede definir como un kernel de transición [12], creado por una serie de distribuciones condicionales en un esquema markoviano, actualizado entre ellas mismas; estas se basan en la distribución a posteriori £(X), donde A es el parámetro para producir una cadena de Markov cíclica alrededor de las condicionales y que converja a dicha distribución.

El siguiente es un resumen del proceso iterativo del muestreador de Gibbs, si únicamente se requiere hacerlo para el parámetro A:

Elegir valor inicial: X^[0]

En la j ésima iteración el ciclo se completa generando valores de cada parámetro condicionados a los valores previamente encontrados de las k distribuciones dadas por X^[]~§[X\X^lj-1]].

Así, se incrementa j hasta la convergencia de la distribución a posteriori del parámetro [12]. Cuando existen más parámetros, cada uno se actualiza a partir del cambio de los anteriores, por ejemplo:

3. METODOLOGíA

En este artículo se estudian cuatro métodos estadísticos: tres modelos ba-yesianos y uno de generación de variable aleatoria. En las tres subsecciones siguientes se describe el proceso analítico de cada una de las tres técnicas: método de inferencia predictiva bayesiana, método de valor esperado baye-siano de la función predictiva y valor esperado bayesiano con parámetros autorregresivos, y en la siguiente, el algoritmo Tabú para optimizar las últimas dos; procesos necesarios para mostrar la metodología estadística que lleva a cada ecuación usada en los pronósticos, cuyo desempeño específico se muestra en los resultados para el caso de estudio al que fueron aplicados. Estas se comparan con otra técnica clásica: la generación de variable aleatoria Poisson, por medio del MAPE, el cual muestra mejor desempeño cuando su valor es menor.

El caso de estudio se trata de un producto de alto volumen de ventas en la empresa, que constantemente busca mejores técnicas para optimizar sus pronósticos, y con ello, la planeación de sus inventarios. La empresa proporciona 24 datos de demanda mensual real. A continuación se explica el proceso analítico de los métodos bayesianos propuestos.

Pronósticos con la función predictiva bayesiana

Los datos tienen una naturaleza discreta que de acuerdo con una prueba de bondad de ajuste siguen una distribución de Poisson, cuyo parámetro X se asume con distribución a priori Gamma, que a su vez tiene parámetros a y fi que se asumirán con distribución a priori uniforme. La distribución Gamma es recomendada en [8]. La función predictiva será explicada a continuación.

La media de los datos con distribución Poisson es 100.17, que se establece como valor inicial para la función de verosimilitud dada por L(X|datos), donde X representa una tasa de ventas mensuales de naturaleza continua. A partir de esta se obtuvo la distribución a posteriori §(X | a, fi, datos) condicional, con la cual se realiza la función predictiva, como se expresa de forma analítica a continuación.

Función de Verosimilitud: Distribución Poisson para los datos

Función de distribución a priori para Lambda: dicho parámetro seguirá la distribución Gamma (a, /?). Para simplificar las ecuaciones se usará la palabra: Datos para representar la muestra: X₁, X₂,...X_n.

Para las ecuaciones (5) y (6), el parámetro X se actualiza usando su estimador, que es el valor esperado de la distribución a posteriori Gammas+£ _x,, p + n). A partir de la distribución dada en (4), actualizándola con los parámetros vía MCMC por las distribuciones dadas por (5) y (6), se estima la función predictiva para pronosticar X_t+1:

Tomando , y agregando términos en la integral y externamente para completar la forma de la distribución gamma cuya integral vale 1, la expresión anterior se puede simplificar así:

A partir de la función predictiva dada en (10) se simulará el valor futuro X ₊₁, o pronóstico de las ventas para este caso, usando los siguientes pasos:

Pronóstico con el valor esperado

El valor esperado para la predicción bayesiana se deduce a partir de la a posteriori obtenida y es condicional a los datos pasados, así E[X_t+₁ |x], donde X_t+1 es el valor que se va a pronosticar en el periodo expos (t+1) y

187

x representa los datos pasados, lo que por definición se ilustra de manera analítica a continuación:

Distribución a priori para X: Gamma (a, fifi), haciendo un cambio en la forma de escritura: 1/fi en lugar de fi.

La distribución a posteriori para X, asumiendo los parámetros a y fi fijos, es:

La integral de la ecuación (13) corresponde a la forma de una distribución Gamma, como se hizo en la sección anterior, agregando términos tanto en la integral como externamente; la expresión (13) se puede simplificar así:

La ecuación (16) tiene la forma de una distribución . El valor esperado para la predicción bayesiana se basa en la distribución a posteriori para X, como se muestra a continuación:

Valor esperado:

En la ecuación (20) la integral tiene la forma de distribución Gamma; agregando términos de forma similar a la expresión (16) se tiene:

El valor esperado de la ecuación (21) constituye la base para elaborar el pronóstico de la demanda de cada mes. Con base en esta se realiza un método de optimización Tabú, que ha mostrado tener óptimos resultados en [17], [19]. Para esta técnica, los parámetros a_ty fí_t se convierten en variables de decisión y se optimizan 24 valores, 12 para cada uno, correspondientes a los 12 meses del año 2011.

Pronóstico con parámetros autorregresivos del valor esperado

Los pronósticos incorporan valores actualizados hasta t-1 de los datos; por lo tanto, en esta sección se propone incorporar una modificación al valor esperado (21), agregando retardos de sus parámetros, los cuales se adicionan multiplicándolos por el coeficiente 0.1, el cual sale del análisis de autocorrelación de orden 1 que tienen los datos originales, por tanto, se toma como un efecto que puede causar dicho retardo en la estimación siguiente.

La ecuación (22) es el valor esperado modificado con dichos retardos:

Donde

X+₁: valor que se va a pronosticar en el tiempo t+1.

n : cantidad de datos disponibles para pronosticar.

a_f: parámetro de forma para el período t y a_t1 es su retardo de orden 1.

f)_t: parámetro de escala para el periodo t y ¡3_t1 es su retardo de orden 1.

A partir de la ecuación (22), para el valor esperado bayesiano se aplica la técnica Tabú para encontrar los parámetros respectivos que minimizan el MAPE.

Metaheurística Tabú

Este es un método que busca optimizar el valor de una función objetivo f(x) con n variables de decisión, de manera que se genere una lista con los valores de los puntos en el espacio n+1 dimensional (n variables y el valor de la función objetivo f(x)), llamada lista Tabú. Dicha lista representa la memoria del recorrido de la búsqueda, para evitar repetir los puntos evaluados previamente, hasta encontrar el mejor valor entre los puntos explorados. En este trabajo se propone una variación al presentado en [20], quienes se basan en una explotación de puntos dentro de un vecindario y una exploración para buscar la mejor solución dentro del espacio de posibles soluciones, que no garantiza encontrar el óptimo global.

La variación propuesta en este trabajo consiste en la exploración de puntos nuevos usando una combinación de una distribución uniforme y una normal, que no se propone en [20], de manera que se encuentren parámetros mensuales para minimizar el error de pronóstico de un producto terminado, aplicado a un caso empresarial. La lista Tabú se actualiza cada que se encuentra un punto que genera un valor de la función objetivo menor al mínimo de todos los demás que están en la lista Tabú.

Dicho algoritmo se aplica al método del valor esperado bayesiano de las ecuaciones (21) y (22), buscando minimizar el MAPE (ecuación 23).

Donde e_t=X_t -X_m es el error entre el valor real (X_t) y el pronóstico (X_t+1), es decir, la función objetivo es el promedio del error relativo absoluto de 12 pronósticos con respecto a los valores reales respectivos de 2011.

Variables de decisión: 24 parámetros: 12 para a_ty 12 para f_t (t=1,2,.....12).

El proceso de estimación usando el valor esperado en cualquiera de los dos casos mostrados en las ecuaciones (21) o (22) es el siguiente:

- Estimar un valor esperado inicial con la sumatoria y valor de n de los primeros 12 meses del año 2010 y asignar un vector inicial a_t y f_t.

- Actualizar la demanda del siguiente año (2011), cambiando la suma-toria y el valor de n en cada periodo de tiempo t o mes.

- Realizar la búsqueda estadística del vecindario de valores de la función objetivo, variando los 24 parámetros a y f , como variables de decisión, y almacenar en la lista del algoritmo Tabú los mejores puntos que minimizan el MAPE y posteriormente la depuran, tratando de encontrar un valor muy cercano al óptimo

- Estimar cada error relativo y promediar los de todo el año (12) para obtener el MAPE de pronósticos de 2011.

4. RESULTADOSY DISCUSIóN

En esta sección se detalla la prueba de bondad de ajuste para los datos, la generación de variable aleatoria y la aplicación de las ecuaciones y métodos bayesianos para encontrar los 12 pronósticos del año 2011, así como la medición del indicador MAPE respectivo a cada método, comparando los resultados para el caso de estudio.

Se estimó una prueba de bondad de ajuste a 12 valores de las ventas mensuales reales del producto elegido, del año 2010 (nombre reservado por confidencialidad). Al realizar la prueba Chi cuadrado, el Valor-P resultante es 0.3154>0.05; lo anterior prueba que los datos se ajustan a una distribución Poisson con el parámetro lambda estimado X. De acuerdo con este análisis, la distribución tiene una media de 100,17 y una desviación de 10,01.

Pronóstico con variable aleatoria

La distribución Poisson con X=100.17 es una alternativa de pronóstico cuando no se cuenta con muchos datos históricos, mientras un modelo de regresión o un ARIMA no lo serían debido a la exigencia de una cantidad grande de datos.

Con fines de ilustración, al usar la distribución Poisson, luego de 1000 simulaciones, se encuentra un error relativo absoluto medio (MAPE para 12 meses) del 23.3 % con un rango de variación entre 18.7 y 27.5 %, como se muestra en la tabla 1, lo cual representa un alto error para efectos de una planeación de producción adecuada para este producto.

Pronósticos con la función predictiva bayesiana

El proceso analítico de la técnica de inferencia predictiva bayesiana explicado previamente se aplica a los datos usando una programación en el lenguaje R, usando la ecuación (10). Proceso que finaliza con el pronóstico de 12 meses de demanda del producto para el año 2011, con los cuales se calculan los errores relativos absolutos, llevando al promedio: MAPE. Este indicador fue simulado 1000 veces, y se encontró un rango de variación entre 27,9 y 58,9 %.

La simulación basada en la ecuación (10) y el muestreador de Gibbs mostró que los pronósticos encontraron incorrelación en los valores simulados, lo cual muestra que hay convergencia; a pesar de esto y de acuerdo con su MAPE, no genera óptimos resultados para el caso de estudio. Se elaboró una prueba de bondad de ajuste en relación con el parámetro simulado X, y se encontró la forma de la distribución Gamma, como originalmente se propuso, lo cual mostró eficiencia del algoritmo de muestreo diseñado.

Simulación de escenarios

Se simularon 12 valores de ventas, a partir de tres escenarios: procesos AR1 con correlaciones 0.1, 0.5 y 0.8, con la media de 100.17. El proceso se repitió 500 veces, y los resultados son muy similares al caso de estudio, acorde con la tabla 3.

La media de errores absolutos relativos MAPE es muy similar en los tres casos.

Pronóstico del Valor esperado con optimización Tabú

En esta sección se utiliza el valor esperado bayesiano ilustrado en la ecuación (22). Es de anotar que sin emplear la metaheurística de optimización Tabú, el indicador MAPE de pronóstico oscila entre 17 y 20 %, pero cuando se utiliza los resultados mejoran significativamente, reduciendo el MAPE a 13.21 %, como se observa en la tabla 4, que además ilustra los resultados obtenidos para los parámetros a, ¡3 que mejor ajustan el pronóstico mensual.

Simulación de escenarios

Se utilizaron los tres escenarios del caso anterior. Los resultados se presentan en la tabla 5. Se aprecia la mejora en el error medio al utilizar el método de optimización Tabú.

En la siguiente subsección se muestra que al incorporar los retardos de cada parámetro el MAPE mejora aun más.

Parámetro autorregresivo del valor esperado con optimización Tabú

Al incorporar el retardo de orden 1 de cada parámetro, ecuación (22), el MAPE se reduce a 8,86 %, lo cual muestra una mejora considerable en los pronósticos con respecto a la metodología anterior.

Aunque este sea una aplicación específica, sí muestra que el algoritmo de búsqueda metaheurístico Tabú mejora mucho el desempeño del método bayesiano por sí solo, aun cuando se incorporan dependencias pasadas, lo cual da ventajas para el uso de esta técnica de pronósticos propuesta.

La tabla 6 ilustra los mejores resultados de los parámetros a y ¡ para cada tiempo encontrados con el algoritmo Tabú diseñado e incorporando el retardado de orden 1 en los parámetros, para el caso real.

Simulación de escenarios

Utilizando los mismos escenarios descritos previamente, los resultados de MAPE se aprecian en la tabla 7. De todos los resultados mostrados para pronosticar, estos son los mejores.

Síntesis de resultados

El indicador que mide el error porcentual medio en el pronóstico (MAPE) para cada uno de los modelos aplicados se muestra en la tabla 8, tanto para el caso real de ventas del producto como para los casos simulados.

En la tabla 8 claramente se observa que el mejor resultado, tanto para el caso real estudiado como para los escenarios de simulación, es el valor esperado con retardos en los parámetros y uso de optimización Tabú, por obtener el menor error medio MAPE de 12 pronósticos.

5. DISCUSIóN

Es de anotar que el método bayesiano debe actualizarse en cada tiempo t, así como el valor de la distribución a priori, de manera que se renueva su pronóstico a partir de información reciente; esto, sumado al uso del algoritmo Tabú, muestra para este caso un mejor acierto en el pronóstico de este producto terminado.

Extensiones de este trabajo pueden llevar a buscar híbridos de heurísticas de optimización en busca de un mejor acierto en la búsqueda de los parámetros, así como permite proponer otros modelos, empleando otras distribuciones, como la normal para el parámetro de la Poisson, o usar métodos de elicitación o, incluso, modelos de regresión bayesiana, o modelos lineales dinámicos en este campo.

CONCLUSIONES

Cuando existen pocos datos no puede usarse una técnica estadística clásica como un modelo ARIMA; si bien para algunos casos se usan solo distribuciones de probabilidad para generar variable aleatoria, esto no considera dependencias a corto o largo plazo o los impactos de una variación temporal, que lo hacen poco acertado, como lo muestra este caso de estudio.

Las técnicas propuestas en este trabajo alrededor de la teoría bayesiana proporcionan alternativas para pronosticar a corto plazo; pero al ser combinadas con el algoritmo Tabú muestran mejor desempeño, tanto para el caso real sobre las ventas del producto como para los escenarios simulados. Lo anterior sugiere que sí es posible proponer diferentes técnicas de pronóstico con relación a los modelos clásicos cuando existen pocos datos, con una precisión adecuada.

La metaheurística Tabú usada es una alternativa de optimización adecuada para una función objetivo que requiera alta precisión. Sin embargo, esta no garantiza un óptimo global.

Por otro lado, la programación de las técnicas bayesianas presentadas en este trabajo no requiere sistemas de cómputo muy complejos o largos periodos de simulación; estas pueden hacerse con el software libre R.

REFERENCIAS

[I] D. Simchi-Levi, P. Kaminski, and E. Simchi-Levi, Designing and Managing the Supply Chair, 3^rd ed., 2008, p. 498.

[2] B. Sani and B. Kingsman, "Selecting the best periodic inventory control and demand forecasting methods for low demand items", J. Oper. Res., vol. 48, n° 7, pp. 700-713, 1997.

[3] V. Gutiérrez and C. Vidal, "Modelos de Gestión de Inventarios en Cadenas de Abastecimiento: Revisión de la Literatura", Rev. Fac. Ing. la Univ., pp. 134149, 2008.

[4] H. Sarimveis, P. Patrinos, C. D. Tarantilis, and C. T. Kiranoudis, "Dynamic modeling and control of supply chain systems: A review", Comput. Oper. Res., vol. 35, n° 11, pp. 3530-3561, Nov. 2008.

[5] R. Watson, "The effects of demand-forecast fluctuations on customer service and inventory cost when demand is lumpy", J. Oper. Res. Soc., vol. 38, n° 1, pp. 75-82, 1987.

[6] F. Diebold, Elementos de pronósticos. México: International Thomson editores, 1999.

[7] R. Cohen and F. Dunford, "Forecasting for Inventory Control: An Example of When 'Simple' Means 'Better'", Interfaces (Providence), vol. 16, n° 6, pp. 9599, 1986.

[8] R. Neelamegham and P. Chintagunta, "A Bayesian model to forecast new product performance in domestic and international markets", Mark. Sci., vol. 18, n° 2, pp. 115-136, 1999.

[9] C. Pedroza, "A Bayesian forecasting model: predicting U.S. male mortality", Biostatistics, vol. 7, n° 4, pp. 530-550, Oct. 2006.

[10] Q. Flora Lu, "Bayesian Forecasting of Stock Prices Via the Ohlson Model," WORCESTER POLYTECHNIC INSTITUTE, 2005. Thesis-Degree of Master of Science in Applied Statistics. Retrieved from https://www.wpi.edu/Pubs/ETD/Available/etd-050605-155155/unrestricted/Flora_Thesis_May_2005.pdf

[II] P. M. Yelland, "Bayesian forecasting of parts demand", Int. J. Forecast., vol. 26, n° 2, pp. 374-396, Apr. 2010.

[12] J. Gill, Bayesian Methods-A social and Behavioral Sciences Approach, 2002, p. 459.

[13] C. J. Barrera and J. C. Correa, "Distribución predictiva bayesiana para modelos de pruebas de vida vía MCMC"/"The Bayesian Predictive Distribution in Life Testing Models via", Rev. Colomb. Estadística, vol. 31, n° 2, pp. 145-155, 2008.

[14] E. Flores, "Teoría de Credibilidad en el Cálculo de Reservas", in Memorias Coloquio e Estadística, 2012.

[15] J. Harrison and C. Stevens, "Bayesian Forecasting", J. R. Stat. Soc., vol. 38, n° 3, pp. 205-247, 1976.

[16] A. Martin, K. Quinn, and J. H. Park, "MCMC pack: Markov Chain Monte Carlo in R", J. Stat. Softw, vol. 42, n° 9, pp. 1-21, 2011.

[17] A. Urrea and F. Torres, "Optimización de una política de inventarios por medio de búsqueda Tabú", in III Congreso colombiano y I Conferencia Andina internacional, 2006, p. 8.

[18] E. A. Silver, "An overview of heuristic solution methods", J. Oper. Res. Soc., vol. 55, n° 9, pp. 936-956, May 2004.

[19] J. Van den Bergh, J. Belién, P. De Bruecker, E. Demeulemeester, and L. De Boeck, "Personnel scheduling: A literature review", Eur. J. Oper. Res., vol. 226, n° 3, pp. 367-385, May 2013.

[20] R. Chelouah and P. Siarry, "A hybrid method combining continuous tabu search and Nelder-Mead simplex algorithms for the global optimization of multiminima functions", Eur. J. Oper. Res., vol. 161, n° 3, pp. 636-654, March 2005.