📘 Introducción

Este informe analiza el comportamiento del sector agrícola colombiano, enfocándose en la industria azucarera, una actividad clave en regiones como el Valle del Cauca. Para ello, se usa como ejemplo la empresa ficticia Ingenio Azucarero del Valle S.A.S., que depende directamente de la producción de caña, su procesamiento y las condiciones del mercado exterior.

Se seleccionaron tres variables mensuales entre enero de 2012 y diciembre de 2024 para estudiar la dinámica del sector:

  1. Producción de azúcar: muestra la cantidad de azúcar procesada por la industria.
  2. Molienda de caña: representa la materia prima utilizada (caña de azúcar).
  3. Exportaciones de azúcar: indica cuánto del producto se vende al exterior.

Estas variables permiten entender la operación, el desempeño comercial y la relación entre producción e internacionalización del sector.

#Cargar librerías necesarias
library(readxl)  # Para leer archivos Excel
library(tseries)  # Para pruebas de estacionariedad
library(forecast)  # Para modelado ARIMA y pronósticos
library(ggplot2)  # Para visualización de datos
library(plotly)  # Para gráficos interactivos
library(timetk)   #timetk simplifica y acelera el análisis exploratorio, visualización, y preparación de datos temporales para modelado. Es ideal para quienes trabajan con series temporales en un flujo de trabajo "tidy" y buscan integrar análisis visuales, detección de patrones y forecasting en un solo paquete.

Cargar base de datos

library(readxl)
data_col <- read_excel("~/Documents/TALLER 2/Base Caso2.xlsx", col_types = c("date", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric", "numeric", "numeric", "numeric", 
    "numeric"))

⚙️ Metodología

El análisis se realizó mediante herramientas estadísticas de series temporales, usando principalmente tres enfoques:

  1. Descomposición STL (Seasonal-Trend-Loess)
    Se usó para dividir cada serie en tres componentes:

    • Tendencia: comportamiento general en el tiempo.
    • Estacionalidad: patrones que se repiten cada año.
    • Residuo: variaciones aleatorias o no explicadas.

  2. Tasa de crecimiento anual
    Se calculó la tasa de crecimiento año a año, comparando la serie original con su componente de tendencia. Esto ayuda a detectar periodos de aceleración, desaceleración o crisis.

  3. Modelo SARIMA
    Se aplicó un modelo ARIMA con estacionalidad (SARIMA) para hacer pronósticos sobre la producción de azúcar. Este modelo fue seleccionado automáticamente con base en criterios estadísticos y por su buen ajuste visual y técnico.

PASO INDISPENSABLE: Declarar la (s) variable (s) como serie (s) temporal (es):

Producción de Azúcar

# Convertir/declarar variable 1=AZUCAR en serie de tiempo mensual
variable1_ts <- ts(data_col$AZUCAR, start = c(2012, 1), frequency = 12)

Molienda de Caña de Azúcar

# Convertir/declarar el CAN en serie de tiempo  mensual
variable2_ts <- ts(data_col$CAN, start = c(2012, 1), frequency = 12)

Exportaciones de azúcares y confitería

# Convertir/declarar las exportaciones de azucar y confiteria en serie de tiempo mensual
variable3_ts <- ts(data_col$X_AZU, start = c(2012, 1), frequency = 12)

🔍 Extracción de señales

Gráfico inicial de la Producción de Azúcar

library(ggplot2)
library(plotly)
library(scales)

# Convertir la serie temporal a un vector numérico
data_col$variable1 <- as.numeric(variable1_ts)

# Crear vector de fechas
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col))

# Agregar fechas como columna
data_col$fechas <- fechas

# Crear el gráfico con mejoras visuales
grafico_serie <- ggplot(data_col, aes(x = fechas, y = variable1)) +
  geom_line(color = "grey", linewidth = 0.4) +
  geom_point(color = "black", size = 0.1) +
  scale_y_continuous(labels = label_comma()) +  # Eje Y legible
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +  # Fechas claras en eje X
  ggtitle("Variable 1: Producción de Azúcar") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 10),
    axis.title = element_text(size = 12),
    plot.title = element_text(size = 14, face = "bold")
  )

# Gráfico interactivo
ggplotly(grafico_serie)

Esta serie evidencia que el sector azucarero opera bajo una estructura cíclica, donde los meses de mayor y menor producción tienden a repetirse cada año.

Los cambios abruptos (como los observados en 2020) deben ser investigados a fondo, pues pueden tener impactos significativos en la planeación financiera.

Extracción señales Producción de Azúcar

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var1 <- stl(variable1_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var1 <- data.frame(
  Time = rep(time(variable1_ts), 4),  # Tiempo repetido para cada componente (son 4 componentes)
  Value = c(stl_decomp_var1$time.series[, "seasonal"], 
            stl_decomp_var1$time.series[, "trend"], 
            stl_decomp_var1$time.series[, "remainder"], 
            variable1_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable1_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var1, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la Producción de azúcar",
       x = "Tiempo",
       y = "Valor")

# Convertir a gráfico interactivo con plotly
ggplotly(p)

La serie muestra una fuerte estacionalidad, con picos regulares cada año. La tendencia tiene tres fases claras: crecimiento (2012–2015), estabilidad (2016–2020) y ligera caída hasta 2023, con recuperación en 2024.

Extracción señales Molienda de caña de azúcar

# Cargar librerías necesarias
library(ggplot2)
library(plotly)
library(scales)

# Descomposición de la serie temporal
stl_decomp_var2 <- stl(variable2_ts, s.window = "periodic")

# Crear vector de fechas correspondiente
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var2 <- data.frame(
  Fecha = rep(fechas, 4),  # Fecha repetida para cada componente
  Valor = c(
    stl_decomp_var2$time.series[, "seasonal"], 
    stl_decomp_var2$time.series[, "trend"], 
    stl_decomp_var2$time.series[, "remainder"], 
    as.numeric(variable2_ts)
  ),
  Componente = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable2_ts))
)

# Crear gráfico con mejoras visuales
p <- ggplot(stl_df_var2, aes(x = Fecha, y = Valor, color = Componente)) +
  geom_line(size = 0.6) +
  facet_wrap(~Componente, scales = "free_y", ncol = 1) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +  # Eje X claro
  scale_y_continuous(labels = label_comma()) +  # Eje Y legible (cuando posible)
  theme_minimal() +
  theme(
    strip.text = element_text(size = 11, face = "bold"),
    axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 10),
    axis.title = element_text(size = 12),
    plot.title = element_text(size = 14, face = "bold")
  ) +
  labs(
    title = "Descomposición temporal de la Molienda de caña de azúcar",
    x = "Tiempo",
    y = "Valor"
  )

# Convertir a gráfico interactivo
ggplotly(p)

También presenta una marcada estacionalidad anual. La tendencia muestra crecimiento al inicio, estabilidad posterior y repunte en 2024. La fuerte similitud con la producción confirma la relación directa entre ambas variables.

Extracción señales Exportaciones de azúcares y confiteria

# Cargar librerías necesarias
library(ggplot2)
library(plotly)

# Descomposición de la serie temporal
stl_decomp_var3 <- stl(variable3_ts, s.window = "periodic")

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var3 <- data.frame(
  Time = rep(time(variable3_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var3$time.series[, "seasonal"], 
            stl_decomp_var3$time.series[, "trend"], 
            stl_decomp_var3$time.series[, "remainder"], 
            variable3_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable3_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de las exportaciones de azúcar y confitería",
       x = "Tiempo",
       y = "Valor")

# Convertir a gráfico interactivo con plotly
ggplotly(p)

Hay alta volatilidad en la serie original, con menos estacionalidad visible. La tendencia muestra una caída sostenida hasta 2018, seguida de una recuperación clara desde 2019. La menor regularidad implica que factores externos como comercio internacional influyen más que ciclos agrícolas.

Después de la descomposición temporal de cada variable, se extrae la variable ajustada por estacionalidad para graficarla junto con la serie original:

Se crea la variable1 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]

Se crea la variable2 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]

Se crea la variable3 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]

Ahora si se puede graficar las series originales versus la ajustada por estacionalidad

Gráfico serie original VS ajustada Producción de Azúcar

library(ggplot2)
library(plotly)
library(scales)

# Crear vector de fechas correctamente alineado con la serie
fechas_var1 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))

# Crear data frame para facilitar el mapeo de leyendas
df <- data.frame(
  fechas = fechas_var1,
  original = as.numeric(variable1_ts),
  ajustada = as.numeric(variable1_sa)
)

# Gráfico mejorado con etiquetas legibles
grafico_ajustada_var1 <- ggplot(df) +
  geom_line(aes(x = fechas, y = original, color = "Serie Original"), size = 0.5) +
  geom_line(aes(x = fechas, y = ajustada, color = "Serie Ajustada"), size = 0.6) +
  scale_color_manual(values = c("Serie Original" = "grey", "Serie Ajustada" = "black")) +
  scale_y_continuous(labels = label_comma()) +  # Eje Y legible
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +  # Eje X claro
  ggtitle("Producción de Azúcar: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Toneladas de Producción de Azúcar") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 10),
    axis.title = element_text(size = 12),
    plot.title = element_text(size = 14, face = "bold")
  )

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)

La serie ajustada permite observar claramente la evolución estructural, sin el “ruido” estacional. Se observa una leve caída reciente, pero con recuperación hacia 2024.

Gráfico serie original VS ajustada Molienda de caña de azúcar

# Crear vector de fechas correctamente alineado con la serie
fechas_var2 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var2 <- ggplot() +
  geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Molienda de caña de azúcar:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Tonelads de Molienda de cala de azúcar") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)

La serie ajustada revela que, aunque hay estacionalidad, la molienda se mantiene relativamente constante con una reciente mejora.

Gráfico serie original VS ajustada Exportaciones de azúcares y confiteria

# Crear vector de fechas correctamente alineado con la serie
fechas_var3 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var3 <- ggplot() +
  geom_line(aes(x = fechas_var3, y = variable3_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Exportaciones de azúcares:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Miles de dolares en Expor de azúcares") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)

Hay mucha más irregularidad; la estacionalidad no es tan fuerte. La serie ajustada permite ver una recuperación estructural reciente en exportaciones.

Ahora graficamos serie original vs tendencia

Tendencia Producción de Azúcar

library(ggplot2)
library(plotly)
library(scales)  # <- Necesario para formatear los ejes

# Convertir la serie a un vector numérico
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
  geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  scale_y_continuous(labels = label_comma()) +  # <-- Aquí se mejora el eje Y
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +  # Opcional: años claros en eje X
  ggtitle("Producción de Azúcar: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Toneladas de Producción de azúcar") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
    axis.text.y = element_text(size = 10)
  )

# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var1)

La tendencia muestra una curva clara de crecimiento, estancamiento y luego repunte reciente. Ideal para identificar ciclos de expansión o ajuste productivo.

Tendencia Molienda de caña de azúcar

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Molienda de caña de azúcar: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Tonelads de Molienda de caña de azúcar") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var2)

Similar a producción, con una recuperación clara hacia el final. Permite validar que la capacidad operativa sigue el ritmo productivo.

Tendencia Exportaciones de azúcares y confiteria

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])

# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
  geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Exportaciones de azúcares: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Miles de dolares en Expor de azúcares") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var3)

La tendencia refleja una importante recuperación desde 2019 tras varios años de caída. Es crucial para analizar estrategia de comercio exterior.

Comparación de componentes

Al separar las series en tendencia y estacionalidad, se hizo evidente que:

  • La estacionalidad explica una gran parte del comportamiento de las tres variables.

  • Analizar la serie ajustada por estacionalidad permite ver con mayor claridad el comportamiento real del sector (ya que se elimina la parte que se repite cada año).

Relación entre variables

  • Producción y molienda están fuertemente relacionadas: cuando se muele más caña, se produce más azúcar.
  • Exportaciones también están relacionadas, pero influidas por otros factores como demanda internacional y precios.

Tasa de crecimiento de la serie de tendencia y original para la Producción de Azúcar

*Gráfico variable original y tendencia Producción de Azúcar: tasa de crecimiento anual**

library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Producción de Azúcar: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)
library(ggplot2)
library(plotly)
library(dplyr)

# Crear el data frame
datos_var1 <- data.frame(
  fecha = fechas_corregidas_var1,
  crecimiento = tasa_crecimiento_var1,
  tendencia = tasa_tendencia_var1
)

# Filtrar entre julio de 2022 y diciembre de 2024
datos_filtrados <- datos_var1 %>%
  filter(fecha >= as.Date("2022-07-01") & fecha <= as.Date("2024-12-31"))

# Crear gráfico
grafico_crecimiento_var1 <- ggplot(datos_filtrados) +
  geom_line(aes(x = fecha, y = crecimiento), color = "grey", size = 0.7) +
  geom_line(aes(x = fecha, y = tendencia), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Producción de Azúcar: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

En general, tasa cercana a 0% hasta 2022. Desde 2023, la tendencia muestra crecimiento sostenido, lo cual es alentador para el sector.

Tasa de crecimiento de la serie original vs tendencia: Molienda de caña de azúcar

# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Molienda de caña de azúcar: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)
library(dplyr)
library(ggplot2)
library(plotly)

# Crear el data frame
datos_var2 <- data.frame(
  fecha = fechas_corregidas_var2,
  crecimiento = tasa_crecimiento_var2,
  tendencia = tasa_tendencia_var2
)

# Filtrar entre julio de 2022 y diciembre de 2024
datos_filtrados <- datos_var2 %>%
  filter(fecha >= as.Date("2022-07-01") & fecha <= as.Date("2024-12-31"))

# Crear gráfico
grafico_crecimiento_var2 <- ggplot(datos_filtrados) +
  geom_line(aes(x = fecha, y = crecimiento), color = "grey", size = 0.7) +
  geom_line(aes(x = fecha, y = tendencia), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Molienda de caña de azúcar: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

Comportamiento muy similar a la producción, con leve crecimiento estructural a partir de 2023.

Tasa de crecimiento de la serie original vs tendencia: Exportaciones de azúcares y confiteria

# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var3 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Exportaciones de azúcares: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

Mucho más volátil, con picos y caídas abruptas. La tendencia suavizada muestra mejora progresiva hasta 2023, y una caída leve posterior.

En general, la tasa de crecimiento es positiva en los primeros años, pero muestra caídas en momentos específicos, lo que puede reflejar impactos económicos o climáticos.

  • La producción y molienda tienen tasas muy similares, reafirmando su relación directa.

  • Las exportaciones muestran más altibajos, indicando mayor exposición a condiciones externas.

🤖 Modelo ARIMA

División en conjunto de entrenamiento y prueba para la Producción de Azúcar que es la elegida para pronosticar

Esta variable fue seleccionada para el pronóstico por ser la más representativa del sector.

El código siguiente divide una serie temporal (variable1_ts) en dos subconjuntos:

Conjunto de entrenamiento (train): Datos desde enero de 2012 hasta septiembre de 2024. Conjunto de prueba (test): Datos desde octubre de 2024 hasta diciembre de 2024.

Esto se hace para evaluar el desempeño de modelos de predicción en datos no vistos.

# Esta división idealmente podria se 80%-70% de los datos para entrenamiento y 20%-30% para prueba o test

# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 2024 y  el conjunto de prueba o test: noviembre 2024-diciembre 2024 

train_size <- length(variable1_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable1_ts, end = c(2024, 9))  # Entrenamiento hasta septiembre 2024
test_ts <- window(variable1_ts, start = c(2024, 10))  # Prueba inicia desde oct2024

🧬 Modelo ARIMA automático normal

El modelo que se presenta a continuación es sin tener en cuenta el factor estacional.

Identificación automática del modelo ARIMA

library(forecast)

# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)

# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)
## Series: train_ts 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1        mean
##       1.0236  -0.4696  -0.4409  180678.547
## s.e.  0.1952   0.0973   0.2155    3609.854
## 
## sigma^2 = 1.291e+09:  log likelihood = -1820.2
## AIC=3650.4   AICc=3650.81   BIC=3665.55
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 44.19856 35461.41 28352.91 -13.45831 27.26131 1.407163 -0.01636981

Estimación del modelo identificado automatico y validación de Significancia de coeficientes

library(lmtest)

# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)
## 
## z test of coefficients:
## 
##              Estimate  Std. Error z value  Pr(>|z|)    
## ar1        1.0236e+00  1.9516e-01  5.2451 1.562e-07 ***
## ar2       -4.6957e-01  9.7281e-02 -4.8270 1.386e-06 ***
## ma1       -4.4094e-01  2.1554e-01 -2.0457   0.04078 *  
## intercept  1.8068e+05  3.6099e+03 50.0515 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajuste del modelo ARIMA(2,0,1) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts, 
                order = c(2, 0, 1))  # Especificamos directamente (p=2, d=0, q=1)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1        mean
##       1.0236  -0.4696  -0.4409  180678.547
## s.e.  0.1952   0.0973   0.2155    3609.854
## 
## sigma^2 = 1.291e+09:  log likelihood = -1820.2
## AIC=3650.4   AICc=3650.81   BIC=3665.55
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 44.19856 35461.41 28352.91 -13.45831 27.26131 1.407163 -0.01636981

Validación de residuales o errores del modelo

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_auto)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1) with non-zero mean
## Q* = 173.02, df = 21, p-value < 2.2e-16
## 
## Model df: 3.   Total lags used: 24

Pronóstico modelo ARIMA automático dentro de muestra o en el set de prueba

# Generar pronóstico para el conjunto de prueba
forecast_arima_auto <- forecast(darima_auto, h = length(test_ts))  # Predecir los valores futuros

# Crear dataframe para gráfico interactivo del pronóstico
forecast_data_auto <- data.frame(Tiempo = time(forecast_arima_auto$mean), 
                            Pronostico = as.numeric(forecast_arima_auto$mean),
                            Observado = as.numeric(test_ts))

# Graficar pronóstico junto con los valores observados reales
p4auto <- ggplot(forecast_data_auto, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("variable1")

ggplotly(p4auto)  # Convertir el gráfico en interactivo

Interpretación modelo automatico (2,0,1):El modelo muestra una caída continua en la producción, pero en los datos reales se observa una recuperación. Esto indica que no logra captar el comportamiento estacional de la serie. Por esta razón, se optó por utilizar el modelo SARIMA, que sí considera la estacionalidad y mejora la precisión del pronóstico.

Pronóstico automático dentro del set de prueba como tabla

# Cargar librerías necesarias
library(forecast)
library(dplyr)

# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts))

# Crear un dataframe con los valores observados y pronosticados
forecast_table_auto <- data.frame(
  Tiempo = time(arima_forecast_auto$mean),  # Extraer las fechas del pronóstico
  Observado = as.numeric(test_ts),  # Valores reales
  Pronosticado = as.numeric(arima_forecast_auto$mean)  # Valores pronosticados
)

# Mostrar la tabla
print(forecast_table_auto)
##     Tiempo Observado Pronosticado
## 1 2024.750  211029.5     203222.3
## 2 2024.833  169151.2     182293.7
## 3 2024.917  189902.3     171746.0

Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, es decir enero 2025

Es decir, le sumamos al periodo de prueba una observación más. Es decir, se estan pronosticando 4 observaciones o trimestres.

# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts) + 1)

# Extraer el pronóstico del próximo trimestre
next_month_forecast_auto <- data.frame(
  Tiempo = time(next_forecast_auto$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast_auto$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast_auto)
##     Tiempo Pronostico
## 1 2024.750   203222.3
## 2 2024.833   182293.7
## 3 2024.917   171746.0
## 4 2025.000   170776.3
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2025 = 170776.301662718"

🧮 Modelo SARIMA automático

Este modelo podria ser una solución o mejora al modelo arima tradicional ya que recoge el efecto estacional de las variables, es recomendable por tanto para datos que si tienen un componente estacional fuerte.

El modelo ajustado en este ejemplo es un SARIMA(3,0,2)(1,1,1)[12], lo que significa:

(3,0,2): Parte ARIMA no estacional: 2 términos autorregresivos (AR). 0 diferenciación (d), lo que indica que la serie no fue diferenciada. 1 término de media móvil (MA).

(1,1,1)[12]: Parte estacional con periodicidad 12 (mensual si los datos son mensuales): 1 término autorregresivo estacional (SAR). 0 diferenciaciones estacionales. 0 términos de media móvil estacionales (SMA).

El modelo SARIMA(3,0,2)(1,1,1)[12] sugiere que:

  • El ajuste es adecuado según los criterios AIC y BIC.

Identificación automática del modelo SARIMA

# Identificación automática modelo SARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)
## Series: train_ts 
## ARIMA(3,0,2)(1,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2     sar1     sma1      drift
##       1.2492  -0.5393  0.0609  -0.9689  0.5619  -0.1805  -0.8113  -120.5336
## s.e.  0.3221   0.2869  0.1633   0.3016  0.1710   0.1060   0.1057    84.7997
## 
## sigma^2 = 337485123:  log likelihood = -1588.89
## AIC=3195.79   AICc=3197.16   BIC=3222.32

A continuación, se crea el objeto darima para luego poder graficar los valores reales y observados:

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(3,0,2)(1,1,1)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(3, 0, 2),  # (p,d,q) -> (3,0,2)
                seasonal = list(order = c(1, 1, 1),  # (P,D,Q) -> (1,1,1)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(3,0,2)(1,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2     sar1     sma1
##       1.1998  -0.4946  0.0816  -0.9148  0.5334  -0.1885  -0.7787
## s.e.  0.3047   0.2977  0.1655   0.2868  0.1883   0.1060   0.0962
## 
## sigma^2 = 342601106:  log likelihood = -1589.72
## AIC=3195.45   AICc=3196.54   BIC=3219.04
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set -649.8222 17322.12 12184.05 -6.115283 12.32966 0.604698
##                     ACF1
## Training set 0.004093373

Validación de residuales del modelo automatico SARIMA

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,0,2)(1,1,1)[12]
## Q* = 26.594, df = 17, p-value = 0.06429
## 
## Model df: 7.   Total lags used: 24

Al revisar los errores del modelo, vemos que están distribuidos de forma equilibrada alrededor de cero. Lo cual nos indica que el modelo funciona bien y que los resultados que genera son confiables.

Pronóstico con el modelo SARIMA dentro del set de prueba-Gráfico líneas

# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts))  # Predecir los valores futuros

# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean), 
                            Pronostico = as.numeric(forecast_arima$mean),
                            Observado = as.numeric(test_ts))

# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("Unidad Variable 1")

ggplotly(p4)  # Convertir el gráfico en interactivo

El modelo sí anticipa correctamente el comportamiento en “V” de la producción al cierre de 2024. Las diferencias entre lo estimado y lo real son pequeñas, lo que valida la elección del modelo SARIMA.

El modelo funcionó correctamente. Al comparar los valores reales y los estimados en el conjunto de prueba, las dos líneas del gráfico están muy cerca, lo que indica que el modelo sí logró capturar bien el comportamiento de la producción.

También se revisaron los errores del modelo y no mostraron patrones raros, lo cual es buena señal.

Pronóstico del modelo automático SARIMA en el set de prueba-Tabla

# Cargar librerías necesarias
library(forecast)
library(dplyr)

# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(auto_arima_model, h = length(test_ts))

# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
  Tiempo = time(arima_forecast$mean),  # Extraer las fechas del pronóstico
  Observado = as.numeric(test_ts),  # Valores reales
  Pronosticado = as.numeric(arima_forecast$mean)  # Valores pronosticados
)

# Mostrar la tabla
print(forecast_table)
##     Tiempo Observado Pronosticado
## 1 2024.750  211029.5     201666.3
## 2 2024.833  169151.2     144767.3
## 3 2024.917  189902.3     184203.1

📈 Pronóstico Enero 2025

El modelo proyecta que en enero de 2025, la producción de azúcar estará dentro del rango normal para ese mes, siguiendo el mismo patrón estacional que se ha observado en años anteriores.

# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)

# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
  Tiempo = time(next_forecast$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast)
##     Tiempo Pronostico
## 1 2024.750   201666.3
## 2 2024.833   144767.3
## 3 2024.917   184203.1
## 4 2025.000   166162.7
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2025 = 166162.742652186"

💡 Recomendaciones

Con base en el análisis, estas son algunas acciones que podría tomar la empresa Ingenio Azucarero del Valle S.A.S.:

  1. Planificación de recursos
    Usar el pronóstico para organizar la producción y preparar el personal y maquinaria antes de los picos de actividad.

  2. Aprovechar el ciclo exportador
    Planear las exportaciones en los meses donde hay más producción y mejores condiciones de mercado.

  3. Prepararse ante cambios económicos
    La tasa de crecimiento permite detectar señales de alerta (como desaceleraciones). La empresa puede anticiparse y tomar decisiones preventivas.

  4. Incluir pronósticos en la gestión
    Los modelos como SARIMA deberían integrarse a los procesos de planeación financiera, logística y comercial de forma continua.

✅ Conclusión

La industria azucarera colombiana tiene una dinámica estacional muy clara, con relaciones estrechas entre producción, molienda y exportaciones. El análisis con STL, tasa de crecimiento y modelo SARIMA permite entender mejor cómo funciona el sector y anticipar comportamientos futuros. Esto ofrece ventajas clave para una empresa que quiera planear de forma estratégica y reducir incertidumbre.