Universidad Nacional de Colombia
📧 jesdiaz@unal.edu.co
Asignatura: Series de tiempo
Fecha: 25 de mayo de 2026
El presente trabajo tiene como objetivo analizar el comportamiento temporal del balance azucarero en Colombia a partir de una base de datos extraída de la plataforma Agronet, la cual contiene información mensual correspondiente al periodo comprendido entre los años 2000 y 2023. Para el estudio se consideraron dos variables principales: la producción total de azúcar y las ventas de los ingenios al mercado nacional, construyéndose dos series de tiempo independientes con el fin de evaluar su evolución, comportamiento y dinámica temporal.
Inicialmente, se realizó un análisis descriptivo de ambas series mediante medidas estadísticas y representaciones gráficas que permitieron identificar fluctuaciones importantes a lo largo del tiempo. Posteriormente, las variables fueron transformadas en series temporales mensuales utilizando el software estadístico R, facilitando el análisis de sus componentes de tendencia, variabilidad y estacionalidad.
Los resultados evidenciaron comportamientos estacionales asociados a ciclos productivos y de consumo dentro del sector azucarero colombiano. Además, se aplicaron pruebas de estacionariedad y modelos SARIMA para modelar la dependencia temporal de las series y describir su comportamiento dinámico. Los modelos seleccionados presentaron un ajuste adecuado según los criterios estadísticos utilizados y las pruebas diagnósticas aplicadas.
Finalmente, el estudio demuestra la utilidad de las técnicas de series de tiempo implementadas en R para el análisis del balance azucarero, proporcionando herramientas relevantes para comprender la evolución de la producción y comercialización del azúcar en Colombia, así como para realizar proyecciones de corto plazo que apoyen la toma de decisiones en el sector.
Series temporales, modelo SARIMA, producción azucarera, ventas nacionales, estacionariedad, diferenciación, pronóstico, autocorrelación, estacionalidad, Box-Ljung, prueba ADF, prueba KPSS, residuos, sector agroindustrial, análisis estadístico, modelación temporal, Colombia, azúcar.
La agroindustria de la caña de azúcar constituye uno de los sectores agrícolas y económicos más importantes de Colombia, debido a su aporte en la producción nacional, generación de empleo y participación en el comercio internacional. Actualmente, el país produce anualmente cerca de 2.2 millones de toneladas de azúcar, destinando aproximadamente el 60% al abastecimiento del mercado nacional y entre el 30% y 40% a la exportación hacia diferentes países de América y otros mercados internacionales. Además, esta actividad representa una fuente importante de ingresos y desarrollo regional, especialmente en el valle geográfico del río Cauca, donde gran parte de la población depende directa e indirectamente de la producción azucarera.
El sector azucarero colombiano también desempeña un papel relevante en la economía nacional, aportando aproximadamente el 0.6% del Producto Interno Bruto (PIB) y cerca del 3.6% del PIB agrícola. Asimismo, genera alrededor de 286.000 empleos directos e indirectos y beneficia a más de 1.2 millones de personas, consolidándose como una de las principales actividades agroindustriales del país. Adicionalmente, la industria participa activamente en la producción de bioetanol utilizado para la oxigenación de combustibles, ampliando así su impacto dentro de la economía y el sector energético nacional.
En este contexto, el presente trabajo tiene como finalidad analizar el comportamiento temporal de la producción azucarera en Colombia a partir de información mensual obtenida de la plataforma Agronet correspondiente al periodo comprendido entre los años 2000 y 2023. Para ello, se construyeron dos series de tiempo asociadas a la producción total de azúcar y a las ventas de los ingenios al mercado nacional, con el propósito de estudiar sus patrones de comportamiento, tendencia, variabilidad y estacionalidad mediante técnicas de análisis de series temporales implementadas en R.
Finalmente, mediante la aplicación de pruebas de estacionariedad, análisis gráficos y modelos SARIMA, se busca describir la dinámica temporal de las series y generar pronósticos de corto plazo que permitan comprender la evolución de la producción y comercialización del azúcar en Colombia, aportando herramientas útiles para el análisis y la toma de decisiones dentro del sector agroindustrial
La información utilizada en este estudio fue obtenida de la plataforma Agronet, la cual proporciona estadísticas oficiales del sector agropecuario colombiano. Se trabajó con una base de datos mensual correspondiente al periodo comprendido entre los años 2000 y 2023, relacionada con la producción azucarera en Colombia.
A partir de dicha información se construyó una serie temporal asociada a la producción total de azúcar, con frecuencia mensual. Posteriormente, los datos fueron organizados y procesados mediante el software estadístico R para realizar el análisis exploratorio y el modelamiento estadístico de la serie.
Inicialmente, se efectuó un análisis descriptivo mediante gráficos y medidas estadísticas con el fin de identificar patrones de comportamiento, tendencia y estacionalidad. Luego, se evaluó la estacionariedad de la serie utilizando pruebas estadísticas y análisis gráficos de autocorrelación (ACF) y autocorrelación parcial (PACF). Debido a que la serie original no presentaba estacionariedad, se aplicó una diferenciación de primer orden para estabilizar su media.
Posteriormente, se ajustó un modelo SARIMA para representar la dinámica temporal de la serie. La selección del modelo se realizó considerando los criterios de información y el comportamiento de las funciones de autocorrelación. Finalmente, la validación del modelo se llevó a cabo mediante el análisis de residuos y la prueba de Box-Ljung, verificando ausencia de autocorrelación significativa y comportamiento aproximado de ruido blanco en los errores.
Finalmente, con el modelo seleccionado se realizaron pronósticos para los siguientes 12 meses, permitiendo analizar el comportamiento esperado de la producción azucarera en el corto plazo.
El análisis de la serie temporal de producción azucarera se realizó con el propósito de identificar patrones de comportamiento, dependencia temporal y componentes estacionales presentes en los datos mensuales comprendidos entre los años 2000 y 2023. Inicialmente, se efectuó una representación gráfica de la serie original para evaluar visualmente la existencia de tendencia, variabilidad y posibles comportamientos estacionales.
Posteriormente, se aplicaron herramientas de análisis exploratorio mediante las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF), las cuales permitieron estudiar la dependencia temporal de la serie y orientar la identificación del modelo estadístico apropiado. Además, se realizaron pruebas de estacionariedad con el fin de verificar si la serie mantenía media y varianza constantes en el tiempo.
Debido a que la serie original no presentó estacionariedad, se aplicó una diferenciación de primer orden para estabilizar su comportamiento. Luego de esta transformación, se observó una mejora en las propiedades estadísticas de la serie, permitiendo proceder con el ajuste del modelo SARIMA.
El modelo seleccionado fue un SARIMA(1,0,1)(0,1,1)[12] el cual permitió representar adecuadamente la dinámica temporal y la componente estacional mensual de la producción azucarera. La validación del modelo se realizó mediante el análisis de residuos, histogramas, funciones ACF residuales y la prueba de Box-Ljung, verificando ausencia de autocorrelación significativa y comportamiento aproximado de ruido blanco en los errores.
Finalmente, se generaron pronósticos para los siguientes 12 meses con el fin de analizar el comportamiento esperado de la producción azucarera y evaluar la incertidumbre asociada a las predicciones futuras mediante intervalos de confianza al 80% y 95%
Existen cuatro componentes de las series temporales que son la tendencia, la estacionalidad, los ciclos irregulares y la variación. El análisis de series temporales puede aislar cada componente y cuantificar el grado en que cada componente influye en la forma de los datos observados y la previsión puede proyectar el patrón subyacente hacia el futuro. (Bee Dagum E, Bianconcini S.(2016)). Así mismos Los gráficos de series temporales pueden revelar patrones como aleatorios, tendencias, cambios de nivel, períodos o ciclos, observaciones inusuales o una combinación de patrones.(Montgomery DC, Jennings C, Kulahci M. (2015)). El componente de tendencia de una serie temporal muestra la dirección general a largo plazo de los datos, ya sea una tendencia descendente o ascendente en cada período, de manera predecible. La tendencia de los datos puede ser lineal o no lineal, dependiendo de las variables consideradas. El componente estacional existe cuando la serie presenta fluctuaciones regulares basadas en las estaciones. La variación estacional ocurre en un período específico, como mensual, trimestral o anual. El componente cíclico, también conocido como ciclos de irregularidad, muestra movimientos oscilatorios de la tendencia en la serie temporal que ocurren durante más de un año. Los datos muestran aumentos y disminuciones que no son periódicos y se repiten a lo largo de un período de tiempo prolongado. Cualquier variación que no se explique por los componentes anteriores (tendencia, estacionalidad y cíclico) se denomina componente aleatorio o irregular. Las perturbaciones no son predecibles, ya que tienden a no seguir la tendencia general exhibida en los datos de la serie temporal.
La estacionariedad es una propiedad fundamental que debe satisfacer una serie de tiempo para poder ser modelada mediante la familia de modelos ARIMA y SARIMA. Una serie se considera estacionaria cuando su media, varianza y estructura de covarianza permanecen constantes a lo largo del tiempo, lo que implica que la serie no presenta tendencia sistemática ni cambios en su variabilidad. Cuando una serie no cumple esta condición, como ocurrió con la serie original analizada, es necesario aplicar transformaciones como la diferenciación para inducir la estacionariedad y garantizar la validez de los modelos estimados.
El supuesto de normalidad establece que los residuos del modelo estimado deben seguir una distribución normal con media cero y varianza constante, expresado como \(ε_t∼N(0,σ^2)\). Este supuesto es relevante para garantizar la validez de las inferencias estadísticas sobre los parámetros estimados, como los intervalos de confianza y las pruebas de significancia. Para verificarlo se aplican pruebas formales como Jarque-Bera o Shapiro-Wilk sobre los residuos del modelo, donde un valor p superior a 0.05 indica que no hay evidencia suficiente para rechazar la hipótesis de normalidad. Generalmente, los histogramas, diagramas de tallo y hojas, diagramas de caja, gráficos de porcentaje-porcentaje gráficos de cuantil-cuantil,gráficos de la función de distribución acumulativa empírica y otras variantes de gráficos de probabilidad son los más útiles para verificar el supuesto de normalidad.
La independencia, también conocida como ausencia de autocorrelación en los residuos, establece que los errores del modelo no deben estar relacionados entre sí. Es decir, el error cometido en un período no debe influir ni predecir el error del período siguiente. Cuando este supuesto se cumple, se dice que los residuos se comportan como ruido blanco, lo que indica que el modelo ha capturado adecuadamente toda la estructura dinámica de la serie y no queda información sistemática sin explicar. Para verificar este supuesto se aplica la prueba de Ljung-Box, cuya hipótesis nula establece que no existe autocorrelación en los residuos hasta un determinado número de rezagos,además el correlograma ACF de los residuos permite una verificación visual, donde la ausencia de barras significativas fuera de las bandas de confianza corrobora que los residuos son independientes y el modelo está correctamente especificado.
La varianza de los residuos debe ser constante, lo cual se puede comprobar mediante un diagrama de dispersión. Este diagrama debe presentar una forma rectangular alrededor del nivel horizontal cero, sin mostrar ninguna tendencia.
\(\textbf{Modelo Autorregresivo de orden p AR(p)}\)
El modelo autorregresivo de orden \(p\), denotado AR(p), generaliza el modelo AR(1) al incluir \(p\) rezagos pasados de la serie para explicar su valor actual. Este modelo es adecuado cuando la serie presenta dependencia con múltiples períodos anteriores y el correlograma PACF muestra un corte abrupto en el rezago \(p\). Su estructura permite capturar dinámicas más complejas de persistencia en la serie, siendo especialmente útil en series económicas y financieras donde los efectos de períodos pasados se prolongan varios períodos hacia adelante.
\(X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_p X_{t-p} + \varepsilon_t\)
\(X_t - \phi_1 X_{t-1} - \phi_2 X_{t-2} - \cdots - \phi_p X_{t-p} = \varepsilon_t\)
\(X_t - \phi_1 B X_t - \phi_2 B^2 X_t - \cdots - \phi_p B^p X_t = \varepsilon_t\)
\((1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p) X_t = \varepsilon_t\)
\(\textbf{Modelo de Media Móvil de orden q MA(q)}\)
El modelo de media móvil de orden \(q\), denotado MA(q), generaliza el MA(1) al considerar los efectos de \(q\) perturbaciones aleatorias pasadas sobre el valor actual de la serie. Este modelo es adecuado cuando la serie presenta memoria de corto plazo en sus errores y el correlograma ACF muestra un corte abrupto en el rezago \(q\). Su capacidad para capturar el efecto acumulado de múltiples choques aleatorios lo convierte en una herramienta valiosa para modelar series donde las perturbaciones externas tienen un impacto que se prolonga varios períodos antes de desvanecerse.
\(X_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}\)
\(X_t = (\varepsilon_t + \theta_1 B \varepsilon_t + \cdots + \theta_q B^q \varepsilon_t)\)
\(X_t = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t\)
\(\textbf{No estacionariedad}\)
Cuando una serie de tiempo no es estacionaria en media, es necesario aplicar el operador de diferenciación \((1-B)^d\) para transformarla en una serie estacionaria. Este operador calcula las diferencias sucesivas entre observaciones consecutivas, eliminando la tendencia presente en la serie. El orden de diferenciación \(d\) indica cuántas veces debe diferenciarse la serie para alcanzar la estacionariedad, siendo \(d=1\) el caso más común en la práctica. De manera análoga, el operador \((1-B^s)^D\) realiza la diferenciación estacional, eliminando los patrones repetitivos asociados al ciclo \(s\) de la serie.
\((1-B)^d\)
\((1 - B^s)^D X_t\)
\(\textbf{modelo Modelo Autorregresivo Estacional AR(P)}\)
El modelo autorregresivo estacional AR(P) es una extensión del modelo AR clásico que captura dependencias en los rezagos estacionales de la serie. En este modelo, el valor actual de la serie se explica mediante sus propios valores pasados en períodos separados por el ciclo estacional \(s\), lo que permite modelar patrones que se repiten regularmente en el tiempo, como comportamientos mensuales o trimestrales. Se utiliza cuando la serie presenta una estructura de dependencia entre observaciones del mismoperíodo en años o ciclos anteriores.
\[X_t = \Phi_1 X_{t-s} + \Phi_2 X_{t-2s} + \cdots + \Phi_P X_{t-Ps} + \varepsilon_t\]
\[X_t - \Phi_1 X_{t-s} - \Phi_2 X_{t-2s} - \cdots - \Phi_P X_{t-Ps} = \varepsilon_t\]
\[X_t - \Phi_1 B^s X_t - \cdots - \Phi_P B^{Ps} X_t = \varepsilon_t\]
\[(1 - \Phi_1 B^s - \cdots - \Phi_P B^{Ps}) X_t = \varepsilon_t\]
\(\textbf{Modelo de Media Móvil Estacional MA(Q)}\)
El modelo de media móvil estacional MA(Q) modela el valor actual de la serie como una combinación lineal de los errores o choques aleatorios ocurridos en períodos estacionales anteriores. A diferencia del AR, no utiliza valores pasados de la propia serie sino los errores pasados en rezagos múltiplos del período estacional \(s\). Es especialmente útil para capturar el efecto de perturbaciones aleatorias que tienen un impacto prolongado a lo largo de los ciclos de la serie.
\[X_t = \varepsilon_t + \Theta_1 \varepsilon_{t-s} + \cdots + \Theta_Q \varepsilon_{t-Qs}\]
\[X_t = (1 - \Theta_1 B^s + \cdots + \Theta_Q B^{Qs})\varepsilon_t\]
\[\Downarrow\] \(\textbf{Modelo SARIMA completo}\)
El modelo SARIMA (Seasonal Autoregressive Integrated Moving Average) es el modelo más general de esta familia, ya que integra en una sola estructura los componentes autorregresivos, de media móvil y de diferenciación, tanto en su parte regular como en su parte estacional. Este modelo es ampliamente utilizado en el análisis y pronóstico de series de tiempo que presentan simultaneamente tendencia, estacionalidad y comportamiento estocástico, permitiendo capturar toda la dinámica de la serie mediante un único modelo parsimonioso. Su notación general es \(SARIMA(p,d,q)(P,D,Q)[s]\), donde los parámetros minúsculos corresponden a la parte no estacional y los mayúsculos a la parte estacional.
\[(1 - \phi_1 B - \cdots - \phi_p B^p)(1 - \Phi_1 B^s - \cdots - \Phi_P B^{Ps})(1-B)^d(1-B^s)^D X_t\] \[= (1 + \theta_1 B + \cdots + \theta_q B^q)(1 + \Theta_1 B^s + \cdots + \Theta_Q B^{Qs})\varepsilon_t\]
\[\Downarrow\]
\[\Phi(B)\, \phi_P(B)\, (1-B)^d\,
(1-B^s)^D\, X_t = \Theta(B)^s\, \theta(B)\, \varepsilon_t\]
La estimación de parámetros postula que existen varios métodos, como el de momentos, el de máxima verosimilitud y el de mínimos cuadrados, que pueden emplearse para estimar los parámetros del modelo identificado tentativamente. por otra parte, es el proceso mediante el cual se calculan los valores numéricos de los coeficientes del modelo SARIMA a partir de los datos observados de la serie. En otras palabras, es el paso donde el modelo deja de ser teórico y se convierte en una ecuación con números concretos.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11041 165419 192019 187140 217716 263434
Desviación estándar
## [1] 38445.57
El análisis descriptivo de la variable Producción Total revela que los datos presentan una producción promedio de 187,140 toneladas, con una mediana de 192,019 toneladas, siendo esta última ligeramente superior a la media, lo que sugiere una leve asimetría negativa en la distribución causada por la presencia de valores extremadamente bajos en la serie. La desviación estándar de 38,445 toneladas, equivalente aproximadamente al 20.5% de la media, refleja una variabilidad moderada en los niveles de producción a lo largo del período analizado. El 50% central de las observaciones se concentra entre 165,419 y 217,716 toneladas, según lo indica el rango intercuartílico. Destaca el valor mínimo de 11,041 toneladas, considerablemente alejado del primer cuartil, lo que confirma la presencia de al menos un valor atípico extremo en la serie, consistente con la caída abrupta observada en el gráfico de la serie histórica, mientras que el valor máximo de 263,434 toneladas se encuentra dentro de un rango coherente con el comportamiento general de la producción.
El histograma evidencia que la mayor parte de la producción se concentra entre 150000 y 240000 toneladas durante el período analizado. La distribución presenta una ligera asimetría hacia la izquierda. Sin embargo, la mayoría de las observaciones se agrupa alrededor de la media, reflejando una variabilidad moderada en la serie
El boxplot de la producción total muestra que la mayoría de los datos se concentran entre aproximadamente 165000 y 220000 toneladas, con una mediana cercana a las 190000 toneladas. Además, se observan algunos valores atípicos inferiores, indicando meses donde la producción fue considerablemente más baja de lo habitual. La variabilidad de la producción es moderada y la mayor parte de los datos se concentra cerca de la mediana.
##
## Augmented Dickey-Fuller Test
##
## data: serie_1
## Dickey-Fuller = -7.4472, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
##
## KPSS Test for Level Stationarity
##
## data: serie_1
## KPSS Level = 0.48506, Truncation lag parameter = 5, p-value = 0.04503
Para evaluar la estacionariedad de la serie temporal de producción total se aplicaron las pruebas ADF y KPSS, debido a que ambas permiten analizar la presencia de raíz unitaria y verificar la estabilidad de la serie.
La prueba ADF tiene como hipótesis.
\(H_{0}:\)La serie de tiempo tiene una raíz unitaria
\(H_{1}:\) La serie es estacionaria.
Arrojo un p-valor menor al nivel de significancia, indicando evidencia de estacionariedad.sin embargo, la prueba KPSS, con hipótesis
\(H_{0}:\) La serie temporal es estacionaria
\(H_{1}:\) La serie no es estacionaria,
presentó un valor-p de 0.04 menos al nivel de significancia, sugiriendo que la serie aún conserva componentes no estacionarios.
lo cual es una señal de comportamiento límite de la serie. Esto implica la posible presencia de componentes estructurales como tendencia o estacionalidad leve.
por esta razón se aplico una diferencia a la serie.
la serie diferenciada oscila de forma relativamente estable alrededor de cero, lo que confirma que la diferenciación fue efectiva para eliminar la tendencia y que la serie es estacionaria, aunque los picos extremos que se observan alrededor de 2009-2010, corresponden probablemente a eventos atípicos externos en la serie, como por ejemplo crisis financiera.
##
## Augmented Dickey-Fuller Test
##
## data: serie_dif
## Dickey-Fuller = -5.6897, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
##
## KPSS Test for Level Stationarity
##
## data: serie_dif
## KPSS Level = 0.044505, Truncation lag parameter = 5, p-value = 0.1
podemos observar que los resultados de los test, indican que la serie diferenciada es estacionaria. La prueba ADF rechaza la hipótesis nula de no estacionariedad (p-value = 0.01), mientras que la prueba KPSS no rechaza la hipótesis nula de estacionariedad (p-value = 0.1). En conjunto, ambas pruebas confirman que la serie puede considerarse estacionaria,
Con base en los correlogramas de la serie diferenciada, el ACF (Figura 5) muestra autocorrelaciones positivas y significativas en los primeros rezagos con un decaimiento gradual,indicando la presencia de un componente autorregresivo, además de un pico negativo pronunciado y significativo en el rezago estacional (aproximadamente −0.4), indicando un término de media móvil estacional de orden 1 (SMA(1)). Por su parte, el PACF (Figura 6) presenta un pico positivo altamente significativo en el primer rezago (≈0.45) seguido de un corte abrupto, confirmando un componente AR(1) en la parte no estacional, mientras que el pico significativo observado en el rezago estacional apunta a un componente autorregresivo estacional de orden 1 (SAR(1)).
Figura 7 descomposición de la serie
La descomposición aditiva evidencia claramente los componentes de tendencia, estacionalidad y ruido.. La tendencia revela un comportamiento no lineal con una fase de crecimiento entre 2000 y 2007, una caída pronunciada que podría asociarse a choques externos, y fluctuaciones posteriores con tendencia decreciente al final del período. La componente estacional muestra un patrón cíclico estable lo que confirma la presencia de estacionalidad, y los residuos se comportan predominantemente como ruido aleatorio, aunque se identifican algunos valores atípicos en los años aproximadamente 2009 y 2023.
## Series: serie_1
## ARIMA(1,0,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ma1 sma1
## 0.7708 -0.3540 -0.7627
## s.e. 0.0712 0.1068 0.0408
##
## sigma^2 = 4.6e+08: log likelihood = -3102.58
## AIC=6213.17 AICc=6213.32 BIC=6227.59
A partir del análisis de los correlogramas ACF y PACF de la serie diferenciada, y con base en las pruebas de estacionariedad que confirmaron la necesidad de una diferenciación estacional de orden 1, se procedió a identificar el modelo mediante la función auto.arima aplicada a la serie original. El modelo seleccionado fue el ARIMA(1,0,1)(0,1,1)[12], que en notación SARIMA indica: un componente autorregresivo de orden 1 (p=1) y un componente de media móvil de orden 1 (q=1) en la parte no estacional, sin necesidad de diferenciación regular (d=0) dado que la serie es estacionaria en media; y en la parte estacional, una diferenciación de orden 1 (D=1) junto con un componente de media móvil estacional de orden 1 (Q=1), con un periodo s=12 que confirma la estacionalidad mensual de la serie. La selección de este modelo se justifica por presentar los menores valores de AIC = 6213.17 y BIC = 6227.59.
\(\textbf{Como nuestro modelo es un SARIMA (1,0,1)(0,1,1)[12] nos queda que:}\)
\(\textbf{operador de rezago}\)
\(B X_t = X_{t-1}\)
\(\textbf{Modelo Autorregresivo de orden 1 AR(1)}\)
\[X_t = \phi_1 X_{t-1} + \varepsilon_t\]
\[X_t - \phi_1 X_{t-1} = \varepsilon_t\]
\[X_t - \phi_1 B X_t = \varepsilon_t\]
\[(1 - \phi_1 B) X_t = \varepsilon_t\]
\(\textbf{Modelo de Media Móvil de orden 1 MA(1)}\)
\[X_t = \varepsilon_t - \theta_1 \varepsilon_{t-1}\]
\[X_t = (1 - \theta_1 B)\varepsilon_t\]
\(\textbf{Diferenciación estacional}\)
\((1-B^s)^D\) con \(s=12\) y \(D=1\),
\[(1 - B^{12})^1 X_t = X_t - X_{t-12}\] \(\textbf{Media móvil estacional de orden 1 SMA(1)}\)
\[X_t = \varepsilon_t - \Theta_1 \varepsilon_{t-12}\]
\[X_t = (1 - \Theta_1 B^{12})\varepsilon_t\]
\(\textbf{Modelo SARIMA(1,0,1)(0,1,1)[12] completo}\)
\[(1 - \phi_1 B)(1 - B^{12}) X_t = (1 - \theta_1 B)(1 - \Theta_1 B^{12})\varepsilon_t\]
## Series: serie_1
## ARIMA(1,0,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ma1 sma1
## 0.7708 -0.3540 -0.7627
## s.e. 0.0712 0.1068 0.0408
##
## sigma^2 = 4.6e+08: log likelihood = -3102.58
## AIC=6213.17 AICc=6213.32 BIC=6227.59
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -885.1456 20872.9 14410.56 -5.163154 12.06564 0.6515802 0.01103693
Mediante la función auto.arima, se identificó el modelo ARIMA(1,0,1)(0,1,1)[12] como el más adecuado según los criterios de información AIC (6213.17) y BIC (6227.59). La parte no estacional del modelo incluye un componente autorregresivo de orden 1 (ar1 = 0.7708), que indica una fuerte dependencia entre el valor actual y el período inmediatamente anterior, y un componente de media móvil de orden 1 (ma1 = −0.354), que corrige parcialmente el error del período previo. En la parte estacional, el modelo incorpora una diferenciación estacional de orden 1 y un término de media móvil estacional (sma1 = −0.7627), reflejando una fuerte corrección de los errores estacionales del año anterior. Todos los coeficientes resultaron estadísticamente significativos. En cuanto al ajuste, el modelo presenta un MAPE de 12.07%, indicando un error porcentual moderado, y un ACF1 de 0.011, muy cercano a cero, lo que sugiere que los residuos no presentan autocorrelación significativa y se comportan aproximadamente como ruido blanco.
La gráfica muestra los residuos del modelo a través del tiempo y se usa para verificar si los residuos son aproximadamente estacionarios.
La gráfica de los residuos muestra que estos fluctúan alrededor de cero y no presentan una tendencia sistemática en el tiempo,indicando que los residuos mantienen un comportamiento estable alrededor de cero. Además, la variabilidad parece mantenerse relativamente constante, aunque se observan algunos valores atípicos extremos. En general, no se evidencian patrones claros de autocorrelación, por lo que el modelo parece capturar adecuadamente la estructura de la serie temporal.
##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.9288, p-value = 2.021e-10
El p-value = 6.635e-06 es extremadamente pequeño, por lo tanto se rechaza la hipótesis nula de normalidad. Los residuos del modelo no siguen una distribución normal, lo cual puede estar influenciado por los valores atípicos extremos que se observaron en la serie.
El histograma muestra una distribución aproximadamente simétrica, aunque con valores extremos en ambas colas. aunque la distribución de los residuos es aproximadamente simétrica y centrada en cero, presenta colas pesadas hacia ambos extremos, lo que evidencia la presencia de valores atípicos que alejan la distribución de la normalidad perfecta. La curva de densidad no se ajusta del todo bien a las barras, siendo más estrecha que la distribución real de los residuos
figura 10 grafica de normalidad de los residuos.
El gráfico Q-Q confirma el rechazo de normalidad: aunque los puntos centrales siguen de cerca la línea roja de referencia, las colas se desvían considerablemente, especialmente la cola izquierda donde los valores se alejan drásticamente hacia abajo
Las gráficas ACF y PACF de los residuos del modelo SARIMA(1,0,1)(0,0,1) muestran que la mayoría de las autocorrelaciones se encuentran dentro de las bandas de confianza. Esto indica ausencia de autocorrelación significativa en los residuos, sugiriendo que el modelo ajustado captura adecuadamente la dependencia temporal y estacional de la serie.
La prueba de Box-Ljung se realizó para comprobar la correlación serial de la siguiente manera:
\(H_0\) No hay autocorrelación serial de la serie temporal.
\(H_1\) Existe autocorrelación serial de la serie temporal.
##
## Box-Ljung test
##
## data: residuals(modelo)
## X-squared = 10.728, df = 10, p-value = 0.3791
La prueba de Box-Ljung aplicada a los residuos del modelo arrojó un p-value = 0.3791. Debido a que este valor es mayor que el nivel de significancia de 0.05, no se rechaza la hipótesis nula de independencia de los residuos. Por lo tanto, se concluye que no existe autocorrelación significativa en los residuos y que el modelo ajustado describe adecuadamente la estructura temporal de la serie.
El gráfico de residuos versus valores ajustados permite visualizar el supuesto de homocedasticidad del modelo. Se observa que los residuos se distribuyen predominantemente alrededor de cero a lo largo de toda la escala de valores ajustados,lo que es consistente con el supuesto de media cero de los errores. Sin embargo, se aprecia una ligera reducción en la dispersión de los residuos a medida que aumentan los valores ajustados, con una mayor variabilidad en los valores más bajos de la serie, lo que sugiere la presencia de una leve heterocedasticidad.
Con el fin de estimar el comportamiento futuro de la serie temporal, se realizó un pronóstico utilizando el modelo SARIMA ajustado. Se generaron predicciones para los próximos 12 periodos, con el fin de evaluar la evolución esperada.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 2023 216543.16 189057.76 244028.6 174507.86 258578.5
## Oct 2023 200916.99 171139.37 230694.6 155376.06 246457.9
## Nov 2023 141706.95 110647.46 172766.5 94205.56 189208.4
## Dec 2023 180730.04 148933.39 212526.7 132101.26 229358.8
## Jan 2024 165934.65 133708.02 198161.3 116648.28 215221.0
## Feb 2024 186249.31 153769.90 218728.7 136576.35 235922.3
## Mar 2024 176889.82 144261.16 209518.5 126988.59 226791.1
## Apr 2024 143744.29 111027.27 176461.3 93707.93 193780.7
## May 2024 96320.44 63551.04 129089.8 46203.97 146436.9
## Jun 2024 160396.91 127596.42 193197.4 110232.89 210560.9
## Jul 2024 208235.44 175416.49 241054.4 158043.20 258427.7
## Aug 2024 225582.63 192752.73 258412.5 175373.63 275791.6
El modelo proyecta que la producción mantendrá un comportamiento oscilatorio durante los próximos 12 meses, sin evidenciar una tendencia marcada. Además, los intervalos de confianza aumentan conforme avanza el horizonte temporal, reflejando una mayor incertidumbre en las predicciones futuras.
A partir del análisis de series de tiempo realizado sobre la producción azucarera en Colombia, se identificó inicialmente que la serie original no presentaba estacionariedad, por lo que fue necesario aplicar una diferenciación de primer orden para estabilizar su comportamiento. Posteriormente, mediante el análisis de las funciones ACF y PACF y el ajuste del modelo SARIMA, se logró capturar adecuadamente la dinámica temporal de la serie.
Los diagnósticos realizados sobre los residuos mostraron ausencia de autocorrelación significativa, lo cual fue confirmado mediante la prueba de Box-Ljung, indicando que los residuos se comportan aproximadamente como ruido blanco y que el modelo ajustado es adecuado.
Además, el análisis permitió identificar periodos con fuertes incrementos y disminuciones en la producción azucarera, evidenciando fluctuaciones importantes a lo largo del tiempo. Finalmente, el pronóstico realizado para los siguientes 12 meses sugiere que la producción continuará presentando un comportamiento oscilatorio, sin una tendencia claramente definida, aunque con una incertidumbre creciente conforme aumenta el horizonte de predicción.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 38536 108669 119181 119318 128944 167772
## [1] 15488.87
La variable Ventas_nacional presenta un comportamiento relativamente simétrico, con una media de 119,318 y una mediana de 119,181, lo que indica ausencia de sesgo marcado. La mayor concentración de los datos se encuentra entre 108,669 y 128,944. La desviación estándar de 15,488.87 evidencia una variabilidad moderada alrededor del promedio, mientras que el rango total muestra la existencia de valores extremos en distintos periodos.
El histograma de ventas al mercado nacional (Figura 1) muestra que la producción azucarera colombiana, la distribución exhibe una asimetría negativa generada por una cola izquierda, esta asimetría explica por qué la media resulta ligeramente superior a la mediana en el análisis descriptivo, y sugiere que aunque las ventas suelen mantenerse relativamente estables durante la mayor parte existen períodos con caídas pronunciadas que afectan la distribución general de los datos.
El boxplot de ventas al mercado nacional (Figura 2) muestra que el sector azucarero colombiano comercializa habitualmente entre 110,000 y 130,000 toneladas métricas de azúcar mensualmente dentro del territorio nacional, con una mediana cercana a 120,000 toneladas
La serie de tiempo (Figura 3) presenta una tendencia creciente desde el año 2000 hasta aproximadamente 2015, período tras el cual se estabiliza con fluctuaciones regulares hasta 2021. A lo largo de todo el período la serie mantiene un patrón estacional que se repite regularmente cada año. Hacia el final de la serie se observa una caída abrupta y atípica que se aleja considerablemente del comportamiento habitual. En términos generales la serie no parece estacionaria debido a la tendencia creciente presente en su primera etapa y los cambios en su nivel a lo largo del tiempo, lo que requiere la aplicación de diferenciación para su posterior modelación.
## Warning in adf.test(serie_2): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: serie_2
## Dickey-Fuller = -5.2753, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
## Warning in kpss.test(serie_2): p-value smaller than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: serie_2
## KPSS Level = 1.6554, Truncation lag parameter = 5, p-value = 0.01
Las pruebas formales de estacionariedad aplicadas a la serie corroboraron lo observado visualmente en la Figura 3. La prueba ADF arrojó un valor p = 0.01 rechazando la hipótesis nula de raíz unitaria, mientras que la prueba KPSS obtuvo igualmente un valor p = 0.01 rechazando su hipótesis nula de estacionariedad, generando resultados contradictorios entre ambas pruebas. Ante esta ambigüedad, característica común en series con raíces unitarias débiles o quiebres estructurales como el observado al final del período, se optó por aplicar una diferenciación a la serie como medida conservadora para garantizar su estacionariedad y proceder con la modelación
La serie diferenciada (Figura 4) muestra que tras aplicar una diferenciación de orden 1 la serie oscila alrededor de cero sin presentar tendencia creciente ni decreciente, lo que indica que la media se estabilizó. La varianza se mantiene relativamente constante a lo largo del período,la serie diferenciada presenta un comportamiento estacionario.
## Warning in adf.test(serie_2_diff): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: serie_2_diff
## Dickey-Fuller = -10.38, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
## Warning in kpss.test(serie_2_diff): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: serie_2_diff
## KPSS Level = 0.015029, Truncation lag parameter = 5, p-value = 0.1
Al aplicar las pruebas formales sobre la serie diferenciada ambas convergieron hacia el mismo diagnóstico, confirmando lo observado visualmente en la Figura 4. La prueba ADF rechazó contundentemente la hipótesis nula de raíz unitaria (Dickey-Fuller = −10.38, p = 0.01) y la prueba KPSS no rechazó su hipótesis nula de estacionariedad (KPSS = 0.015, p = 0.10), confirmando que la serie diferenciada es estacionaria. Esto valida que una única diferenciación de orden 1 fue suficiente para inducir la estacionariedad requerida
El ACF muestra un corte rápido después del lag 1, con un pico significativo en el rezago estacional (lag 1.0), indicando un componente de media móvil en la parte no estacional y un componente de media móvil estacional SMA(1) y el PACF presenta múltiples rezagos significativos negativos en la parte inicial con un decaimiento gradual, y un pico significativo en el rezago estacional (lag 1.0), evidenciando un componente autorregresivo en la parte no estacional y un componente autorregresivo estacional SAR.
En conjunto ambos correlogramas sugieren que el modelo candidato para la serie diferenciada de ventas nacionales del sector azucarero es un SARIMA con componentes MA(1) en la parte regular y componentes estacionales en lag 12, cuyos órdenes exactos serán confirmados mediante la función auto.arima.
La descomposición aditiva muestra que la tendencia tiene un crecimiento sostenido entre 2000 y 2015, seguido de fluctuaciones decrecientes hacia el final del período con una caída pronunciada en 2022, la componente estacional presenta un patrón cíclico regular y constante, Los residuos se comportan como ruido aleatorio alrededor de cero.
m<-auto.arima(serie_2)
m
## Series: serie_2
## ARIMA(0,1,3)(0,0,2)[12]
##
## Coefficients:
## ma1 ma2 ma3 sma1 sma2
## -0.6543 -0.1275 -0.1141 0.3825 0.1757
## s.e. 0.0597 0.0726 0.0635 0.0601 0.0554
##
## sigma^2 = 137090140: log likelihood = -3051.8
## AIC=6115.59 AICc=6115.9 BIC=6137.47
Mediante la función auto.arima aplicada a la serie original de ventas nacionales, se identificó el modelo ARIMA(0,1,3)(0,0,2)[12] como el más adecuado. La parte no estacional indica que no fue necesario un componente autorregresivo (p=0), se confirmó la diferenciación de orden 1 (d=1) consistente con el análisis previo, y se incorporaron tres términos de media móvil (q=3). En la parte estacional el modelo no requirió diferenciación ni autorregresión estacional, pero sí incluyó dos términos de media móvil estacional (Q=2) con periodo mensual s=12. El modelo fue seleccionado por presentar los menores valores de AIC = 6115.59 y BIC = 6137.47
\(\textbf{Diferenciación regular (1-B)^1}\)
Dado que la serie de ventas nacionales no era estacionaria en media, fue necesario aplicar el operador de diferenciación regular de orden 1, el cual calcula la diferencia entre el valor actual y el valor del período inmediatamente anterior, eliminando así la tendencia presente en la serie y garantizando la estacionariedad requerida para la modelación.
\[(1-B)^1 X_t = X_t - X_{t-1}\]
\(\textbf{Media Móvil de orden 3 MA(3)}\)
El modelo de media móvil de orden 3, denotado MA(3), expresa el valor actual de la serie como una combinación lineal del error actual y los errores de los tres períodos inmediatamente anteriores. Este modelo es adecuado cuando el correlograma ACF presenta un corte después del tercer rezago, indicando que el impacto de los choques aleatorios se disipa en tres períodos.
\[X_t = \varepsilon_t - \theta_1\varepsilon_{t-1} - \theta_2\varepsilon_{t-2} - \theta_3\varepsilon_{t-3}\]
\[X_t = (1 - \theta_1 B - \theta_2 B^2 - \theta_3 B^3)\varepsilon_t\] \(\textbf{Media móvil estacional de orden 2 MA(2)}\)El modelo de media móvil estacional de orden 2, modela el valor actual de la serie en función de los errores ocurridos en los dos ciclos estacionales anteriores, es decir en el mismo mes de los dos años previos.
\[X_t = \varepsilon_t - \Theta_1\varepsilon_{t-12} - \Theta_2\varepsilon_{t-24}\]
\[X_t = (1 - \Theta_1 B^{12} - \Theta_2 B^{24})\varepsilon_t\]
\(\textbf{Modelo SARIMA(0,1,3)(0,0,2)[12] completo}\)Integrando todos los componentes identificados, el modelo SARIMA(0,1,3)(0,0,2)[12] combina la diferenciación regular de orden 1, el componente de media móvil MA(3) y el componente de media móvil estacional SMA(2), con período estacional \(s=12\) correspondiente a datos mensuales. La expresión general en notación de operadores es:
\[(1-B)(1-B^{12})^0 X_t = (1 - \theta_1 B - \theta_2 B^2 - \theta_3 B^3) (1 - \Theta_1 B^{12} - \Theta_2 B^{24})\varepsilon_t\]
Como \(D=0\) el operador estacional desaparece quedando:
\[(1-B) X_t = (1 - \theta_1 B - \theta_2 B^2 - \theta_3 B^3) (1 - \Theta_1 B^{12} - \Theta_2 B^{24})\varepsilon_t\]
## Series: serie_2
## ARIMA(0,1,3)(0,0,2)[12]
##
## Coefficients:
## ma1 ma2 ma3 sma1 sma2
## -0.6543 -0.1275 -0.1141 0.3825 0.1757
## s.e. 0.0597 0.0726 0.0635 0.0601 0.0554
##
## sigma^2 = 137090140: log likelihood = -3051.8
## AIC=6115.59 AICc=6115.9 BIC=6137.47
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 377.4698 11584.21 8719.566 -0.8115169 7.761446 0.7929609
## ACF1
## Training set -0.0003179086
El modelo ARIMA(0,1,3)(0,0,2)[12] el modelo está compuesto principalmente por términos de media móvil tanto en la parte no estacional como en la estacional. Las estimaciones muestran que los errores recientes tienen mayor influencia que los estacionales. Las métricas de error (MAPE ≈ 7.76%) indican un buen nivel de precisión, y la ausencia de autocorrelación en los residuos (ACF1 ≈ 0) sugiere que el modelo captura adecuadamente la dinámica de la serie.
La Figura 8 muestra que los residuos oscilan alrededor de cero sin tendencia creciente ni decreciente y con varianza relativamente constante a lo largo del período, lo que indica que el modelo capturó adecuadamente la estructura de la serie.
##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.96833, p-value = 6.635e-06
La prueba de Shapiro-Wilk arrojó un valor p = 6.635×10⁻⁶, significativamente menor a 0.05, por lo que se rechaza la hipótesis nula de normalidad, indicando que los residuos no siguen una distribución normal.
El histograma de los residuos (Figura 9) muestra una distribución aproximadamente simétrica y centrada en cero, con una distribución cercana a una forma normal representada por la curva roja. Sin embargo se aprecia una cola izquierda lo cual es consistente con el rechazo de normalidad detectado en la prueba de Shapiro-Wilk.
La Figura 10 muestra que los residuos se distribuyen aleatoriamente alrededor de cero a lo largo de toda la escala de valores ajustados, sin presentar patrones sistemáticos ni tendencia alguna, lo que sugiere un comportamiento homocedástico aceptable.
El gráfico Q-Q muestra que la mayoría de los puntos se alinean sobre la línea roja de referencia en la zona central, lo que indica un comportamiento aproximadamente normal en ese rango. Sin embargo se observan desviaciones en ambas colas, se alejan considerablemente de la línea, confirmando lo detectado en la prueba de Shapiro-Wilk.
Los correlogramas ACF y PACF de los residuos (Figuras 12 y 13) muestran que prácticamente todas las autocorrelaciones se encuentran dentro de las bandas de confianza, lo que indica que los residuos no presentan estructura temporal significativa y se comportan como ruido blanco. Se observan algunos rezagos que rozan levemente los límites de las bandas, especialmente en el PACF, pero sin superarlos de forma contundente
##
## Box-Ljung test
##
## data: residuals(m)
## X-squared = 3.8252, df = 10, p-value = 0.9549
La prueba de Ljung-Box arrojó un valor p = 0.9549, ampliamente superior a 0.05, por lo que no se rechaza la hipótesis nula de independencia, confirmando que los residuos no presentan autocorrelación significativa, este resultado es consistente con lo observado visualmente en los correlogramas ACF y PACF de los residuos.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 2023 122823.5 107818.37 137828.6 99875.14 145771.8
## Oct 2023 122518.1 106641.75 138394.4 98237.32 146798.9
## Nov 2023 118791.6 102581.30 135001.9 94000.10 143583.1
## Dec 2023 116249.3 99964.00 132534.5 91343.10 141155.4
## Jan 2024 117828.3 101468.35 134188.2 92807.94 142848.6
## Feb 2024 114544.8 98110.55 130979.0 89410.80 139678.7
## Mar 2024 116683.9 100175.71 133192.1 91436.81 141931.0
## Apr 2024 114401.8 97819.99 130983.6 89042.11 139761.5
## May 2024 116960.6 100305.44 133615.7 91488.74 142432.4
## Jun 2024 115489.8 98761.69 132217.9 89906.36 141073.3
El gráfico de pronóstico (Figura 14) muestra el comportamiento histórico de las ventas nacionales del sector azucarero colombiano desde el año 2000 hasta 2023, y la proyección generada por el modelo ARIMA(0,1,3)(0,0,2)[12] para los 10 períodos siguientes en color azul. A lo largo del período histórico la serie mantuvo oscilaciones regulares con un patrón estacional mensual definido. La caída atípica extrema registrada hacia 2022 influyó en el punto de partida del pronóstico, razón por la cual el pronóstico inicia desde niveles bajos debido a la caída atípica observada al final de la serie. No obstante el modelo proyecta una recuperación progresiva de las ventas nacionales hacia niveles más consistentes con el comportamiento histórico, con una región sombreada en azul que representa el intervalo de confianza del pronóstico y cuya amplitud se incrementa progresivamente.
El análisis realizado sobre la serie de ventas nacionales permitió identificar un comportamiento estacional marcado y una tendencia creciente durante gran parte del período estudiado. Debido a la falta de estacionariedad observada en la serie original, fue necesario aplicar una diferenciación de primer orden, logrando estabilizar la media y permitiendo posteriormente el ajuste del modelo SARIMA.
El modelo ARIMA(0,1,3)(0,0,2)[12] presentó un ajuste adecuado según los criterios AIC y BIC, además de mostrar residuos sin autocorrelación significativa de acuerdo con la prueba de Ljung-Box. Aunque los residuos no cumplieron completamente el supuesto de normalidad, el modelo logró capturar satisfactoriamente la dinámica temporal y estacional de la serie.
Finalmente, los pronósticos obtenidos sugieren una recuperación gradual de las ventas nacionales luego de la caída atípica observada al final del período analizado, manteniendo un patrón oscilatorio consistente con el comportamiento histórico del sector azucarero colombiano.
De manera general en este trabajo se analizaron dos series temporales relacionadas con el sector azucarero colombiano: la producción total de azúcar y las ventas nacionales. A través del análisis exploratorio, las pruebas de estacionariedad y la modelación SARIMA, fue posible identificar patrones de tendencia, estacionalidad y dependencia temporal presentes en ambas series.
Los resultados mostraron que las dos variables requerían transformaciones para alcanzar estacionariedad, permitiendo posteriormente ajustar modelos capaces de representar adecuadamente su comportamiento dinámico. Los diagnósticos realizados evidenciaron que los modelos seleccionados lograron capturar la mayor parte de la estructura temporal de las series, ya que los residuos se comportaron aproximadamente como ruido blanco y no presentaron autocorrelación significativa.
Además, los pronósticos obtenidos reflejan que tanto la producción como las ventas nacionales continuarían mostrando fluctuaciones estacionales en el corto plazo, aunque con niveles de incertidumbre crecientes conforme aumenta el horizonte de predicción. En términos generales, el análisis confirma que los modelos SARIMA constituyen una herramienta útil para estudiar y proyectar variables del sector agroindustrial colombiano, facilitando la comprensión de su comportamiento y apoyando procesos de planeación y toma de decisiones.
Ministerio de Agricultura y Desarrollo Rural. (2025). Producción de azúcar, alcohol y miel 2000-2023 [Conjunto de datos]. Agronet. https://agronet.gov.co/
Mapuwei, T. W., Ndava, J., Kachaka, M., & Kusotera, B. (2022). Una aplicación del modelo de pronóstico ARIMA de series temporales para predecir la producción de tabaco en Zimbabwe. American Journal of Modeling and Optimization, 9(1), 15-22.
Manoj, K., & Madhu, A. (2014). Una aplicación del modelo de pronóstico ARIMA de series temporales para predecir la producción de caña de azúcar en la India. Facultad de Ciencias Económicas, 9(1), 81-94.
Montgomery, D. C., Jennings, C., & Kulahci, M. (2015). Introducción al análisis y pronóstico de series temporales (2.ª ed.). John Wiley and Sons.
Sector Agroindustrial de la Caña. (s.f.). Asocaña. https://www.asocana.org/
Ruiz-Ramírez, J., Hernández-Rodríguez, G. E., & Zulueta-Rodríguez, R. (2011). Análisis de series de tiempo en el pronóstico de la producción de caña de azúcar. Terra Latinoamericana, 29(1), 103-109. http://www.scielo.org.mx/