Introducción

Hay grandes progresos­ en el desarrollo y las aplicaciones de la predicción del clima a mediano plazo y su predicción estacional (^{Vitart et al. 2012}). Los algoritmos de predicción automáticas más usados son con base en el suavizado exponencial o modelos autorregresivos integrados de media móvil (ARIMA) (^{Hyndman y Khandakar, 2008}). ^{Box y Jenkins (1976)} desarrollaron la metodología clásica que emplea las series de tiempo para generar modelos como el autoregresivo de media móvil (ARMA) o también el modelo ARIMA para obtener predicciones.

^{Karl et al. (2000)} reportaron un aumento en la tasa de calentamiento global usando la serie de tiempo de la temperatura media global indicada por ^{Quayle et al. (1999)}, por medio del análisis de valores mensuales de temperatura y con modelos ARMA. ^{Reikard (2009)} investigó la predicción de la radiación solar en intervalos de tiempos de 5 min hasta varias horas y aunque los datos exhibieron variabilidad no lineal debido a la nubosidad, en la mayoría de las pruebas se obtuvieron los mejores resultados usando los modelos ARIMA. ^{Pulido (2002)} propuso estimar la demanda de agua en las próximas 24 h en un sistema de distribución de agua para riego usando modelos ARIMA y otros modelos. Para predecir la lluvia del monzón de verano en la India, Chatto-padhyay y ^{Chattopadhyay (2010)} identificaron un modelo ARIMA como adecuado, pero el modelo de redes neuronales autorregresivas (ARNN) proporcionó mejores predicciones, mientras que ^{Narayanan et al. (2013)} usaron modelos ARIMA para predecir las lluvias antes del monzón en el oeste de la India.

Debido a que los modelos ARIMA son una herramienta para realizar predicción de series de tiempo univariadas, en esta investigación se propuso elaborar un programa de cómputo que calcule la predicción en tiempo real de variables meteorológicas usando modelos ARIMA y probar su efectividad en condiciones de baja y alta precipitación.

Materiales y Métodos

Para esta investigación se usó una computadora con procesador de 2.2 GHz, 2 GB de memoria RAM y sistema operativo Windows 7®. En la computadora se instaló: el programa de cómputo MySQL Server®, que es un gestor de bases de datos para almacenar información (^{Korhonen et al., 2008}); Microsoft Visual Studio 2010® que es un conjunto completo de herramientas de desarrollo para la generación de aplicaciones Web ASP.NET, Servicios Web XML, aplicaciones de escritorio y aplicaciones móviles (^{Randolph et al., 2010}); MySQL Conector Net 6.3.5® que es un conector del programa Microsoft Visual Studio 2010® con MySQL Server (^{Kofler, 2005}); R Statistics 2.15.3®, un paquete de cómputo estadístico (^{Dalgaard, 2008}); librerías ‘rcom’ y ‘rscproxy’del programa R Statistics 2.15.3 (conectores del programa R Statistics 2.15.3 con Microsoft Visual Studio 2010); y la librería ‘forecast’ del programa R Statistics 2.15.3, se usó para la estimación y predicción de los modelos ARIMA.

Para almacenar información meteorológica se elaboró una base de datos integrada con dos tablas de datos, en el programa MySQL Server (Figura 1). La primera tabla de datos se denominó ‘estación’ y fue usada para guardar la información de cada estación meteorológica, y por cada estación se requiere un identificador de estación, latitud, longitud, altitud y nombre, y la llave primaria es el identificador de estación. La segunda tabla de datos, denominada ‘elemhoraria’, se usó para almacenar la información de los datos meteorológicos a nivel horario de estaciones meteorológicas; los datos almacenados en esta tabla son: fecha y hora, evapotranspiración (ET₀ en mm), velocidad del viento (VELS en m/s), precipitación (mm), radiación solar (RADSOL en W/m²), temperatura media (TEMP en °C), humedad relativa (HR en %), y un identificador de estación de la cual provienen los datos; la llave primaria es la unión de los datos de fecha e identificador de estación.

Figura 1 Base de datos. 

En la Figura 1 se observa que una estación puede tener muchos registros a nivel horario y muchas estaciones pueden tener datos meteorológicos para una hora en particular.

Datos meteorológicos

Para comprobar la bondad predictiva de los modelos ARIMA, se usaron datos de tres estaciones meteorológicas automáticas (EMAS) del Servicio Meteorológico Nacional, México, para el 2013. Las EMAS fueron: ENCB. II del IPN, ubicada en 19° 29’ 55” N, 99° 08’ 43” O y altitud de 2240 m; Acolman, ubicada en 19° 38’ 05” N, 98° 54’ 42” O y altitud de 2269 m; Chapingo, ubicada en 19° 29’ 39” N, 98° 53’ 19” O y altitud de 2260 m.

En las EMAS para este estudio hay datos continuos a nivel horario de cinco variables meteorológicas en dos periodos: el primero es del 7 de marzo de 2013 a las 16:00 h y el 17 de marzo de 2013 a las 15:00 h; el segundo es del 16 de junio de 2013 a las 16:00 h y 26 de junio de 2013 a las 15:00 h. Las variables meteorológicas obtenidas de las EMAS fueron: velocidad del viento (m/s), precipitación (mm), radiación solar (W/m2), temperatura media (°C), humedad relativa (%). Además se calculó la evapotranspiración de referencia (ET0) por el método de Penman Monteith (^{Allen, 2006}) con los datos anteriores.

Modelos ARIMA

Según ^{Pankratz (1983)}, los modelos ARIMA sirven para predecir series simples (de una sola variable), en los que las predicciones de los modelos ARIMA están basadas sólo en valores pasados de la variable a predecir. Los modelos ARIMA se pueden usar para hacer predicciones a corto plazo porque la mayoría de ellos ponen mayor énfasis en el pasado reciente que en el pasado distante; se aplican a variables discretas o continuas, aunque el tiempo debe ser igualmente espaciado y en intervalos discretos; son útiles para predecir series de datos que contienen variación estacional (u otras variaciones periódicas), incluyendo aquellas con patrones estacionales cambiantes; requieren como mínimo alrededor de 50 observaciones; se aplican sólo a series de datos estacionarios, y una serie de tiempo estacionaria tiene una media, varianza y función de autocorrelación constantes a través del tiempo (^{Pankratz, 1983}).

El requisito de una serie de tiempo estacionaria puede parecer enteramente restrictiva, pero la mayoría de las series no estacionarias en la práctica se pueden transformar a una serie estacionaria a través de una transformación llamada “diferenciar”, la cual es una operación relativamente simple que envuelve el cálculo de cambios sucesivos en los valores de las series de datos. Los cambios en la serie de datos se conocen como (w_t) y se obtienen con la ecuación w_t=z_t - z_t-1, donde z representa los valores de la serie de datos. Con las diferencias se construye una nueva serie diferente de la serie original, y una “diferencia” es cuando la media de una serie de datos cambia con el tiempo. Es posible “diferenciar” más de una sola vez para obtener una serie estacionaria. Al ya tener una serie estacionaria, se realiza la búsqueda por un buen modelo ARIMA y consiste en: identificación, estimación, diagnóstico del modelo; y si el modelo es adecuado se realiza la predicción (^{Pankratz, 1983}).

Descripción del procedimiento de la librería Forecast para estimar el modelo ARIMA

Según ^{Hyndman et al. (2013)}, un obstáculo común al usar modelos ARIMA para predecir es que el proceso de selección del orden es generalmente considerado subjetivo y difícil de aplicar. Por tanto, se elaboró la librería Forecast para elegir el orden del modelo de manera automática, y donde los algoritmos son aplicables a ambos, datos estacionales y no-estacionales.

Para ^{Hyndman et al. (2013)} un proceso ARIMA(p,d,q) no-estacional está dado por:

1-Bdyt=c+∅Byt+∅(B)εt

donde {ε_t} es un proceso de ruido blanco con media cero y varianza σ², B es el operador de retraso, y Ø(z) y Ø(z) son polinomios de orden p y q, respectivamente. Para asegurar casualidad e invertibilidad se asume que Ø(z) y Ø(z) no tienen raíces para |z|<1. Si c≠0, hay un polinomio implícito de orden d en la función de predicción. El proceso estacional ARIMA (p,d,q) (P,D,Q)_m está dado por

1-BmD(1-B)dyt=c+∅Bmyt+Θ(Bm)θ(B)εt

donde Ø(z) y Ø(z) son polinomios de orden P y Q respectivamente, cada uno no conteniendo raíces dentro del círculo unitario. Si c≠0, hay un polinomio implícito de orden d+D en la función de predicción.

La tarea principal de la librería Forecast que ^{Hyndman et al. (2013)} realizan en la predicción automática del modelo ARIMA, es seleccionar un apropiado orden de modelo y son los valores de p, q, P, Q, D, d. Si D y d conocidos. Los órdenes p, q, P y Q se pueden seleccionar por medio de un criterio de información como el Criterio de Información de Akaike (AIC):

AIC=-2log⁡L+2(p+q+P+Q+k)

donde k=1 si c≠0 y 0 de otra manera, y L es la probabilidad maximizada del modelo ajustado a los datos diferenciados (1-B^m)^D (1-B)^dy_t.

Para propósitos de predicción ^{Hyndman et al. (2013)} indican que es mejor hacer tan pocas diferencias como sea posible. Para datos no estacionales ^{Hyndman et al. (2013)} consideran modelos ARIMA(p,d,q) donde d es seleccionada basándose en el test de raíces unitarias sucesivas KPSS (^{Kwiatkowski et al., 1992}). El método prueba los datos para una raíz unitaria; si el resultado de la prueba es significativa, se prueban los datos diferenciados para para una raíz unitaria; y así sucesivamente.

Para datos estacionales, en la librería Forecast se consideran modelos ARIMA (p,d,q) (P,D,Q)_m donde m es la frecuencia estacional y D=0 o D=1, dependiendo de una prueba extendida de Canova-Hansen^{(Canova and Hansen, 1995}). Después de seleccionar D se elige d aplicando el test de raíces unitarias sucesivas KPSS a los datos estacionales diferenciados (si D=1) o a los datos originales (si D=0).

Estimación de la predicción

Con Microsoft Visual Studio 2010 se desarrolló una aplicación ejecutable (.exe), con la cual se realizan funciones para la predicción de las variables meteorológicas. Sin embargo, antes de describir dichas funciones es importante remarcar que la mayoría de las EMAS tienen la opción de descargar la información que registran y guardarla en archivos de texto. Por lo anterior, se realizó la primera función para extraer los datos de las variables meteorológicas almacenados en archivos de texto (de cada EMAS) y guardarlos en la base de datos. En la base de datos se almacenan los datos meteorológicos a nivel horario de distintas EMAS y se organizan por fecha y por identificador de EMA; los datos obtenidos se guardan en la tabla de datos ‘elemhoraria’. Con esta función también se calcula la evapotranspiración de referencia por el método de Penman Monteith (^{Allen, 2006}) y el resultado se almacena en la misma tabla de datos.

Cuando termina la primera función se tienen los promedios de las variables meteorológicas a nivel horario en la base de datos. Con la segunda función se genera una serie de tiempo por cada variable meteorológica de las tres EMAS, y la serie de tiempo así generada tiene 60 datos. Cada dato de la serie de tiempo consiste en el promedio de dos horas; por ejemplo, si para el día 16/06/2013 a las 16:00, 17:00, 18:00 y 19:00 h el promedio de temperatura fue 22, 23.7, 24.7 y 24.8 °C respectivamente, y el día 21 de junio de 2013 a las 14:00 y 15:00 h el promedio de temperatura fue 16.2 y 15.6 °C, respectivamente, entonces la serie de tiempo tendrá los valores 22.85, 24.75, …, 15.9 °C, con un total de 60 datos. Las series de tiempo generadas se almacenan en un archivo de texto creado automáticamente con extensión ‘.txt’ para cada EMA (Figura 2); el archivo tiene seis columnas (una columna por cada variable meteorológica) y 61 filas. La primer fila contiene los nombres de las variables meteorológicas; sin embargo, sólo se analizaron cuatro variables. La primera variable está en la primer columna y tiene los datos de evapotranspiración de referencia (mm), en la cuarta columna están los datos de radiación solar (W/m²), en la quinta columna están los datos de la temperatura (°C), y en la sexta columna están los datos de humedad relativa (%).

Figura 2 Contenido del archivo de texto con series de tiempo de variables meteorológicas para la EMA ENC-BII para el periodo de junio. 

Al terminar la segunda función se tienen las series de tiempo para cada variable meteorológica necesaria para realizar la predicción. Con la tercera función se realiza la predicción. En el primer paso de la tercera función se establece una conexión entre el programa R Statistics 2.15.3 con Microsoft Visual Studio 2010; el lenguaje de programación fue C#. Después se envía un comando al programa R Statistics 2.15.3 para hacer el ajuste de las series de tiempo, de cada variable meteorológica, a un modelo ARIMA usando la función “auto.arima” de la librería Forecast (^{Hyndman et al., 2013}); después de lo anterior, ya con el modelo ARIMA estimado automáticamente, se envía un comando al programa R Statistics 2.15 para hacer la predicción de los 60 elementos siguientes en la serie de tiempo.

La función “auto.arima” de la librería Forecast (^{Hyndman et al., 2013}) regresa el mejor modelo ARIMA. No obstante, la función “auto.arima” requiere argumentos como una serie de tiempo univariada, el orden de la primera diferencia “d” (si no se coloca, la función “auto.arima” elige un valor de acuerdo con la prueba KPSS), el orden de la primera diferencia estacional “D” (si no se coloca, la función “auto.arima” lo calcula), el máximo valor para p, q, P, Q, el valor inicial de p, q, P, Q (opcional), si la serie de tiempo estacionaria, si la serie de tiempo es estacional, entre otros datos opcionales.

Para la función “auto.arima” se especificaron las opciones valor máximo de p igual a 5, valor máximo de q igual a 5, y la opción estacional igual a “TRUE” porque de lo contrario se restringe la búsqueda para modelos no-estacionales. La última función consiste en guardar las predicciones obtenidas con la función “auto.arima”. Para lo anterior, se crean archivos de texto con extensión ‘.txt’ y en estos se almacenan las predicciones; el nombre del archivo se crea con la unión del nombre de la EMA, el nombre de la variable y la palabra PRED al final (para indicar que es predicción). En la Figura 3 se muestra el archivo generado de la estación ENCBII para la variable temperatura del aire.

Figura 3 Archivo generado con datos de la predicción estimada con la función “auto.arima” para la EMA ENCBII en el periodo de junio y la variable meteorológica temperatura del aire (°C). 

En el archivo texto generado por la predicción (Figura 3), en la primera línea, la predicción del promedio de la variable entre las 16:00 y 17:00 h el 21 de junio de 2013, en la segunda línea está la predicción del promedio de la variable entre las 18:00 y 19:00 h el 21 de junio de 2013, y en la última fila del archivo está la predicción del promedio de la variable entre las 14:00 y 15:00 h del día 26 de junio de 2013.

Resultados y Discusión

La precisión de la predicción en los dos periodos fue evaluada. En el primero, los datos observados entre 7 de marzo de 2013 a las 16:00 h hasta el 12 de marzo de 2013 a las 15:00 h, se usaron para generar la serie de tiempo, y los datos observados entre 12 de marzo de 2013 a las 16:00 y 17 de marzo de 2013 a las 15:00 h se usaron para compararlos con los datos estimados con el modelo de predicción. Después se evaluó la precisión de la predicción en el segundo periodo. Los datos observados entre el16 de junio de 2013 a las 16:00 h y el 21 de junio de 2013 a las 15:00 h se usaron para generar la serie de tiempo. Los datos observados entre 21 de junio de 2013 a las 16:00 h y el 26 de junio de 2013 a las 15:00 h se usaron para compararlos con los datos de las predicciones estimadas.

En las Figuras 4 y 5 se graficaron los datos observados con los datos estimados de las variables meteorológicas: humedad relativa (%), temperatura del aire (°C), radiación solar (W/m²) y evapotranspiración de referencia (mm). Los datos tomados como observados de evapotranspiración de referencia fueron calculados con el método de Penman Monteith (^{Allen, 2006}) y usando los datos observados de las demás variables meteorológicas.

Figura 4 Predicción y variables observadas 60 h hacia adelante para la EMA Acolman (A), Chapingo (B), y ENCB. II del IPN (C) en el primer periodo (marzo). 

Figura 5 Predicción y variables observadas 60 h hacia adelante para la EMA Acolman (A), Chapingo (B), y ENCB. II del IPN (C) en el segundo periodo (junio). 

Los modelos ARIMA obtenidos para cada estación y variable meteorológica están en el Cuadro 1.

Para una serie de tiempo dada {y_n}, la predicción persistente se obtiene al colocar y(n+1)=y(n), que implica que el promedio de la variable para la siguiente hora es igual al promedio de la variable en la hora actual (^{Kavasseri et al., 2009}).

Para comparar la bondad predictiva de los modelos ARIMA con la predicción persistente se calcularon mediciones del error y del cuadrado medio del error (MSE). ^{Cadenas y Rivera (2007)} mencionan que si el valor observado en el tiempo t es y_t y Ft es la predicción para el mismo tiempo, entonces el error se define como e_t=y_t-F_t, y el cuadrado medio del error es:

MSE=1n∑t=1net2

En nuestro estudio el valor observado del promedio de la temperatura del aire entre las 16:00 y 17:00 del 12 de marzo, 2013, fue 20.7 °C y el valor obtenido con la predicción del modelo ARIMA para el mismo tiempo fue 19.96 °C, por lo cual el MSE del modelo ARIMA, de ese tiempo fue de 0.547. El valor de la predicción persistente de ese mismo tiempo fue 14.95 (el valor observado del promedio de la temperatura del aire entre las 14:00 y 15:00 h del 12 de marzo de 2013); por lo tanto, el valor del MSE del modelo persistente, de ese tiempo fue de 33.06. La misma operación se realizó para las 40 h hacia adelante (20 tiempos hacia delante).

El valor MSE se calculó en tiempos de 5 en 5 hacia adelante para todas las variables (5, 10, 15, 20, etc.) hasta que el valor del MSE obtenido con la predicción de los modelos ARIMA (MSE _A) fuera mayor que el valor del MSE obtenido con la predicción del modelo de persistencia (MSE_B). En la variable temperatura del aire del periodo de marzo de la estación Acolman se encontró que hasta 15 tiempos hacia adelante el valor del MSE_A (2.531) fue menor que el del MSE_B (10.309); y a los 20 tiempos el valor del MSE_A (11.204) fue mayor que el del MSE_B (10.422). Para comparar los errores se calculó un porcentaje de mejoría de la predicción del modelo ARIMA con respecto a la predicción con el modelo persistente.

El porcentaje de mejoría de la predicción del modelo ARIMA con respecto a la predicción con el modelo persistente (PM) se calculó como sigue:

PM=100-MSEA*100MSEp

donde MSE_A es el MSE del modelo ARIMA, MSE_P es el MSE del modelo persistente.

En la variable temperatura del aire del periodo de marzo de la estación Acolman y 15 tiempos hacia adelante, se encontró que el valor del PM fue 75.4 %, lo cual indica que el modelo ARIMA se desempeña 75.4 % mejor que el modelo persistente hasta los 15 tiempos hacia adelante (30 h). De la misma manera se realizó el cálculo del PM para todas las variables meteorológicas en ambos periodos (marzo y junio) y para las tres EMAS (Cuadro 2).

Dos aspectos importantes en un esquema de predicción son: 1) que tan bien un modelo retiene su precisión sobre el horizonte de predicción y, 2) que tan robusto es el esquema para la elección del horizonte de predicción (^{Kavasseri et al., 2009}). Para observar el primer aspecto se realizaron predicciones de los valores de las variables meteorológicas hacia adelante hasta que el PM fuera menor que cero. Cuando los valores de PM son menores que cero, indica que el MSE_P fue mejor que el MSE_A (la predicción con el modelo persistente fue mejor que con el modelo ARIMA). Para observar el segundo aspecto se eligieron dos periodos: el periodo de marzo con muy poca precipitación (menos de tres eventos de precipitación y menos de 1 mm en total), y el periodo de junio en donde se presentan más de 20 eventos de precipitación (más de 8 mm en total).

El modelo ARIMA predice mejor que el persistente más de 15 tiempos hacia adelante para las variables ET0, HR, RADSOL y TEMP en el periodo de marzo (Cuadro 2). Para el periodo de junio se encontró que el modelo persistente fue mejor que el ARIMA para las variables ET0, HR y TEMP. Una probable razón de esto es que la precipitación afecta las demás variables meteorológicas y puede cambiar el comportamiento de una serie temporal.

Conclusiones

El uso de software de computadora y modelos ARIMA permite al investigador estimar la predicción de variables meteorológicas automáticamente y en tiempo real. Sin embargo, los resultados indican que, en promedio, en el periodo de marzo con muy poca precipitación (menos 1 mm), la predicción con los modelos ARIMA fue mejor que la predicción con el modelo persistente en: 53.6 % en evapotranspiración de referencia hasta 20 tiempos hacia adelante (40 h); 46.7 % en humedad relativa hasta 15 tiempos hacia adelante (30 h); 71 % en radiación solar hasta 30 tiempos hacia adelante (60 h); 43.4% en temperatura del aire hasta 15 tiempos hacia adelante (30 h). En el periodo de junio, las predicciones obtenidas con el modelo persistente fueron mejores que con el modelo ARIMA.