INTRODUCCIÓN
El pronóstico de condiciones agrometeorológicas, en sentido general, se refiere a las predicciones esperadas que influyen sobre la producción agropecuaria (Caraza y Quintero, 1997). La utilidad de los mismos se orienta en proveer una guía para la adopción de decisiones, las que en el área de la producción agropecuaria resultan de interés dada la dependencia de esta actividad con las condiciones de la atmósfera.
En la ganadería, se han utilizado las salidas de los modelos globales empleados en la predicción del tiempo asociadas con indicadores como el índice de temperatura y humedad (THI) o rangos de temperatura del aire para la evaluación de las condiciones de los animales (Domínguez et al., 2010). Además de ello, muchas de las técnicas de predicción agrometeorológica se fundan en la relación estadística que existe entre las variables dependientes que han de ser estimadas y las variables agrometeorológicas independientes (precipitación, temperatura), o las variables deducidas (índice de humedad del suelo, influencia de la humedad atmosférica).
La estadística inferencial provee varias herramientas para el análisis que resultan de interés para las ciencias meteorológicas. Los modelos autorregresivos integrados de promedios móviles o ARIMA (acrónimo del inglés autoregressive integrated moving average, también referido como método de Box - Jenkins), donde las estimaciones futuras se sustentan únicamente en los datos previos, han sido aplicados en Cuba para la modelación y pronóstico en diferentes ramas del saber. Por ejemplo, en meteorología se encuentran los trabajos de Osés et al., (2017) y Cabanas (2005), mientras que en la rama agropecuaria se hallan, por ejemplo, los trabajos de Sánchez (2013) y Guevara et al., (2004).
El objetivo general de este trabajo es predecir las condiciones de los animales mediante modelos ARIMA de series temporales para el escenario nacional.
MATERIALES Y MÉTODOS
Se emplearon los datos de 64 estaciones meteorológicas pertenecientes a la red de estaciones de observación del Instituto de Meteorología. Igualmente se usaron los Boletines de Indicadores Fundamentales de la ganadería. De la primera, se utilizaron los datos de la empresa pecuaria La Vitrina.
A continuación se listan los indicadores seleccionados:
Producción de leche
Reproductivos:
Relativos a la masa total:
Para la evaluación del estado de los animales se empleó el Índice de Temperatura y Humedad (ITH) mediante la ecuación 1 (Nienabar y Hahn, 2007; Mujica, 2005 citado por Aguilar, 2006) para un periodo de 21 años (1993 2013).
Donde, t d es la temperatura del bulbo seco (oC) y HR representa la humedad relativa expresada en forma decimal. En la tabla 1 aparecen los intervalos adoptados para la valoración.
Tabla 1.
Escala para la evaluación del ITH.
| Valor del ITH | Evaluación | |
|---|---|---|
| Inferior | Superior | |
| 74 | Normal | |
| 75 | 78 | Alerta |
| 79 | 83 | Peligro |
| 84 | - | Emergencia |
Fuente: Nienabar y Hahn (2007) y Mujica (2005).
Se aplicaron técnicas exploratorias para la identificación de patrones y para ello, los datos fueron divididos y a partir de esto, se crearon series temporales como sigue:
Valores diarios (cada año con 365 datos por estación).
Valores semanales (cada año incluyó 52 datos por estación),
Por valores mensuales (cada año incluyó 12 datos por estación),
Valores anuales (un total de 21 puntos),
Por estaciones (valores diarios, semanales, mensuales y anuales para cada una de las 23 estaciones).
Para el análisis y representación se crearon tres grupos: ITH1, que incluyó las estaciones de las provincias más orientales (desde Camagüey hasta Guantánamo), ITH2 (Matanzas hasta Ciego de Ávila) e ITH3 (Pinar del Rio, Mayabeque, La Habana e Isla de la Juventud).
Se procedió al examen visual de los datos, mediante la representación gráfica con el propósito de identificar irregularidades en las series. Posteriormente se pasó a la descomposición de los mismos (tendencia, ciclicidad y estacionalidad). Se aplicaron los test de Dickey - Fuller y KPSS (Kwiatkowski- Phillips - Schmidt - Shin) para identificar la estacional de las series.
La prueba Dickey - Fuller se basó en asumir que la serie se puede aproximar por un proceso AR(1) con tres variantes: media cero, media diferente de cero y tendencia lineal. En este caso se partió de las hipótesis:
La hipótesis nula (2) asume que la serie no es estacionaria. El estadístico de la prueba se denota por p y su distribución bajo H 0 permite calcular los valores críticos, de tal forma que el criterio de rechazo es p < p0.05, con p es valor calculado del estadístico.
Se representó un diagrama estacional como se describe en Hyndman y Athanasopoulos (2013), es decir, semejante a un gráfico de tiempo excepto que los datos se trazaron contra las estaciones en años separados. La elección de un modelo ARIMA óptimo (p,d,q) se efectuó mediante una función AUTOARIMA. Se consideró a p como el parámetro auto regresivo, d al número de fases de diferenciación no estacionales y q al parámetro de medias móviles.
Los pronósticos de objetos STL (por sus siglas en inglés, Standard Triangle Language) se obtuvieron aplicando un método de predicción no estacional a los datos ajustados estacionalmente y re estadificando utilizando el último año del componente estacional. Se empleó el método Holt Winters de la familia de series temporales de suavizado exponencial.
Para la evaluación e iteración, se plotearon los ACF y PAFC para los residuales de los modelos resultantes. Todos los análisis se efectuaron mediante la aplicación de la herramienta R de R Core Team, (2018).
RESULTADOS Y DISCUSIÓN
Obtención de las series temporales para las diferentes escalas de trabajo
Se adicionaron los datos diarios, semanales y mensuales para conformar series de tiempo de la suma de valores de ITH. Para el ajuste de la serie de tiempo semanal, primeramente se truncaron los primeros y últimos registro de las series y se asignaron los nuevos valores de principio y fin.
Chequeo de la estacionalidad de los datos
Para el ajuste a un modelo tipo ARIMA se requiere que las series sean estacionarias, es decir, que su media, varianza y autocovarianza sean invariantes en el tiempo. En la Tabla 2 aparecen los resultados obtenidos para las pruebas seleccionadas.
Tabla 2.
Test de Dickey - Fuller
Como se puede apreciar en cualquiera de los casos, se rechaza la hipótesis nula y todas las series son estacionarias (p < 0.05). Lo que se corrobora con los resultados del test KPSS (Tabla 3).
Tabla 3.
Test de KPSS para el nivel de estacionalidad.
La autocorrelación en las series es un indicador de cómo cada valor de la serie temporal. En este caso se empleó la función acf() para la detección de esta característica.
El panel en la figura 1 muestra los autocorrelogramas, donde se muestran la correlación entre la serie y sus rezagos. Siempre se tiene que ρ(0)=1. Las líneas discontinuas en azul representan las bandas de confianza de ρ(k) de nivel 95% bajo la hipótesis de que la serie es un ruido blanco (incorrelada).
Figura 1.
Autocorrelogramas para las series.
Se observa que dentro de las series de tiempo hay una fuerte correlación entre cada valor y el valor anterior. Esto significa que los valores de las series temporales no son independientes entre sí. También hay una correlación estadísticamente significativa para el segundo y tercer rezagos. En este caso, el retardo predeterminado es donde N es el número de observaciones y m el número de series.
Una forma simple de lidiar con la autocorrelación es la diferenciación. Diferenciar significa restar valores previos de los valores actuales. La diferenciación de orden 1 significa que se sustrae el valor anterior del valor actual. En este caso, todavía hay autocorrelaciones estadísticamente significativas.
Los gráficos mostrados en la Figura 2 muestran las series temporales.
Figura 2.
Series temporales para los valores medios del ITH.
Una vez obtenidas las series temporales se confeccionaron los diagramas estacionales como se describe en Hyndman y Athanasopoulos (2018).
Los resultados de la aplicación del pronóstico automático mediante suavizado exponencial (método Holt Winters), tanto para la escala semanal como mensual, aparecen en el panel representado en la Figura 3. Para el primer lapso la predicción se realizó para 52 semanas (un año); mientras que en el segundo, se determinó usar 36 meses (tres años).
Figura 3.
Pronóstico automático mediante suavizado exponencial (método Holt Winters)
Para el pronóstico de la serie temporal se empleó la función auto.arima(), mediante la cual se determinaron los valores óptimos (p,d,q) (Hyndman, 2008).
En la Tabla 4 aparecen los estadísticos obtenidos para cada uno de los modelos ARIMA y STL. Los estadísticos Raíz del Error Cuadrático Medio (RMSE), el Error Medio Absoluto (MAE) y el Error Medio Absoluto Porcentual (MAPE), miden la magnitud de los errores, o sea, el mejor será el que muestre el valor más pequeño; y el Error Medio (ME) y el Error Medio Porcentual (MPE) miden el sesgo estadísticamente y se considera que el mejor modelo da un valor próximo a 0 (Sánchez, 2013). Otro criterio lo dan Sánchez y Poveda (2006), para quienes el desempeño de los modelos se evalúa para todas las series mediante el error cuadrático medio (RMS). Sobre los correspondientes en este análisis, tanto para la escala semanal como mensual, se observa que el de mejor ajuste corresponde al grupo oriental: (1,1, 2)(0,0,1)52 y (0,0,1)(0,1,1)12.
Debe considerarse que los modelos ARIMA, por lo general predicen bien para horizontes de tiempo (h) fuera de la muestra cortos y medios. El desempeño de los pronósticos de distintos modelos puede variar según sea la amplitud de dicho horizonte. Cuando se pronostica una serie, el valor de h hasta que sería apropiado pronosticar, estará determinado por la incertidumbre con que se pronostica para cada h (si la amplitud del intervalo de predicción es muy amplia, la información que se brinda será trivial o poco precisa) (Blaconá, 2018).
Tabla 4.
Comparación de modelos probados para la serie índice de temperatura y humedad (ITH).
Para el caso de la aplicación de esta herramienta en el caso particular de una estación meteorológica se escogió la 78308, con una escala temporal diaria para el análisis. En la Figura 4 aparece la serie temporal obtenida al efecto.
Figura 4.
Serie temporal (escala diaria). Estación 78308.
Para la aplicación del pronóstico automático mediante suavizado exponencial (Holt - Winters), como los subsiguientes, requieren el análisis de estacional. La trama estacional de acuerdo con Hyndman y Athanasopoulos (2018) aparece en la Figura 5 a modo de gráfica de tiempo, excepto que los datos se representan en función de las estaciones en años separados.
Figura 5.
Seasonal plot. Estación 78308.
Figura 6.
Autocorrelación para la serie temporal.
Seguidamente se muestran los resultados obtenidos para la predicción automatizada utilizando un suavizado exponencial (método Holt Winters).
Figura 7.
Pronóstico Holt Winters (Estación 78308).
Figura 8.
Pronóstico STL + ETS (A,N,N).
De acuerdo con la literatura consultada, la modelación ARIMA se ha aplicado fundamentalmente en el pronóstico de la producción lechera (Deluyker et al., 1990; Sánchez, 2013). Otros trabajos dentro de esta línea, se aproximan al empleo de una estimación del inventario ganadero bovino, tal y como muestran Cuenca et al. (2008) y Pérez (2014).
CONCLUSIONES
Se logró predecir las condiciones de los animales mediante modelos ARIMA de series temporales para el escenario nacional. Estos métodos de análisis dan una alternativa útil a los modelos lineales clásicos para el análisis de series como las del índice de temperatura y humedad.
Tanto para la escala semanal como mensual, los mejores ajustes se obtuvieron en el grupo correspondiente a la zona oriental de Cuba. Se observa que el mejor corresponde al grupo oriental de Cuba: (1,1, 2)(0,0,1)52 y (0,0,1)(0,1,1)12; con errores medios porcentuales (MPE) igual a -0.017 y -0.06 respectivamente.