#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.
library(lmtest)
library(readxl)
data_col <- read_excel("C:/Users/laura.astudillo/OneDrive - PUJ Cali/Desktop/Maestría en Finanzas/Analitica de negocios - finanzas/Caso #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"))

Introducción

El sector agroindustrial de la caña de azúcar en Colombia es un pilar fundamental de la economía del país, contribuyendo significativamente al PIB y a la generación de empleo. Este sector se destaca por su alta productividad a nivel mundial, con una concentración geográfica en el valle del Cauca. La industria no solo produce azúcar, sino también otros subproductos importantes como etanol y energía.

La industria de la caña de azúcar, pilar fundamental de la economía colombiana, experimentó una notable recuperación en 2024 gracias a condiciones climáticas más favorables. La disminución de las precipitaciones no solo impulsó el rendimiento agrícola, sino que también resultó en un incremento del 6% en la molienda de caña, esto de acuerdo al Informe anual 2024-2025 del Sector Agroindustrial de la Caña (Asocaña), este proceso es un indicador clave de la producción total del sector. Por ello, para el próximo informe financiero, se presenta un análisis detallado de las proyecciones para el mes de enero de 2025, poniendo un énfasis especial en el comportamiento de la molienda de caña, ya que su desempeño histórico y reciente es fundamental para prever la dinámica productiva y comercial del sector.

Metodología

Variable 1 - CAN (Molienda de azúcar)

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

Variable 2 - AZÚCAR (Producción de azúcar)

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

Variable 3 - X_AZU (Exportación de azúcar)

# Convertir/declarar variable 1=ENER 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 variable 1 - CAN (Molienda de azúcar) en niveles -Original

library(ggplot2)
library(plotly)

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

# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)), 
                                      y = variable1)) +
  geom_line(color = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1 - CAN (Molienda de azúcar): Serie original") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  theme_minimal()

ggplotly(grafico_serie)

Extracción señales Variable 1 - CAN (Molienda 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 variable 1 - CAN (Molienda de azúcar)",
       x = "Tiempo",
       y = "Valor")

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

Extracción señales variable 2 - AZÚCAR (Producción de azúcar)

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

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

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

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

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

Extracción señales variable 3 - X_AZU (Exportación de azúcar)

# 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 variable 3 - X_AZU (Exportación de azúcar)",
       x = "Tiempo",
       y = "Valor")

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

Variables ajustadas por estacionalidad

variable 1, CAN (Molienda de azúcar) ajustada por estacionalidad

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

variable 2, AZÚCAR (Producción de azúcar) ajustada por estacionalidad

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

variable 3, X_AZU (Exportación de azúcar) ajustada por estacionalidad

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

Gráficos de las series originales versus la ajustada por estacionalidad

Gráfico serie original VS ajustada Variable 1, CAN (Molienda de azúcar)

# 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))

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var1 <- ggplot() +
  geom_line(aes(x = fechas_var1, y = variable1_ts), color = "#7AC5CD", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 1 - CAN (Molienda de azúcar):
Serie Original vs Serie Ajustada por Estacionalidad")+
  xlab("Tiempo") +
  ylab("Tonelada") +
  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_var1)

Gráfico serie original VS ajustada Variable 2, AZÚCAR (Producción 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 = "#7AC5CD", 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("Variable 2 - AZÚCAR (Producción de azúcar):
Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Tonelada") +
  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)

Gráfico serie original VS ajustada Variable 3, X_AZU (Exportación de azúcar)

# 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 = "#7AC5CD", 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("Variable 3 - X_AZU (Exportación de azúcar):
Serie Original vs Serie Ajustada por Estacionalidad")+
  xlab("Tiempo") +
  ylab("Tonelada") +
  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)

Tendencia Variable 1, CAN (Molienda de azúcar)

library(ggplot2)
library(plotly)

# 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" = "#76EE00", "Tendencia" = "black")) +
  ggtitle("Variable 1 - CAN (Molienda de azúcar): Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Tonelada") +
  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_var1)

Tendencia Variable 2, AZÚCAR (Producción 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" = "#76EE00", "Tendencia" = "black")) +
  ggtitle("Variable 2 - AZÚCAR (Producción de azúcar): 
Serie Original vs Tendencia")+
  xlab("Tiempo") +
  ylab("Tonelada") +
  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)

Tendencia Variable 3, X_AZU (Exportación de azúcar)

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" = "#76EE00", "Tendencia" = "black")) +
  ggtitle("Variable 3 - X_AZU (Exportación de azúcar):
Serie Original vs Tendencia")+
  xlab("Tiempo") +
  ylab("Tonelada") +
  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)

Tasa de crecimiento de la serie de tendencia y original para la variable 1, CAN (Molienda de azúcar)

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var1 <- (variable1_ts[(13:length(variable1_ts))] / variable1_ts[1:(length(variable1_ts) - 12)] - 1) * 100
tasa_tendencia_var1 <- (tendencia_var1[(13:length(tendencia_var1))] / tendencia_var1[1:(length(tendencia_var1) - 12)] - 1) * 100

# Crear vector de fechas corregido, es decir que inicie desde enero 2013
fechas_corregidas_var1 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))

# Verificar longitudes
print(length(fechas_corregidas_var1))
## [1] 144
print(length(tasa_crecimiento_var1))
## [1] 144
print(length(tasa_tendencia_var1))
## [1] 144

Gráfico variable original y tendencia variable 1, CAN (Molienda 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 = "#EEA2AD", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable 1 - CAN (Molienda 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)

# --- Cálculo de la tasa de crecimiento anual ---
tasa_crecimiento_var1 <- (variable1_ts[(13:length(variable1_ts))] / variable1_ts[1:(length(variable1_ts) - 12)] - 1) * 100
tasa_tendencia_var1   <- (tendencia_var1[(13:length(tendencia_var1))] / tendencia_var1[1:(length(tendencia_var1) - 12)] - 1) * 100
# --- Fechas corregidas ---
# La base inicia en 2014-01 → la tasa anual arranca en 2015-01
fechas_corregidas_var1 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))

# --- Filtrar solo enero 2022 - diciembre 2024 ---
inicio <- as.Date("2022-01-01")
fin    <- as.Date("2024-12-01")

filtro <- fechas_corregidas_var1 >= inicio & fechas_corregidas_var1 <= fin

fechas_filtradas      <- fechas_corregidas_var1[filtro]
tasa_crecimiento_filt <- tasa_crecimiento_var1[filtro]
tasa_tendencia_filt   <- tasa_tendencia_var1[filtro]

# --- Gráfico ---
grafico_crecimiento_var1 <- ggplot() +
  geom_line(aes(x = fechas_filtradas, y = tasa_crecimiento_filt), color = "#EEA2AD", size = 0.7) +
  geom_line(aes(x = fechas_filtradas, y = tasa_tendencia_filt), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable 1 - CAN (Molienda de azúcar):
Tasa de crecimiento anual % de la serie Original y la tendencia (2022 - 2024)")+
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

Tasa de crecimiento de la serie de tendencia y original para la variable 2, AZÚCAR (Producción de azúcar)

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))

# Verificar longitudes
print(length(fechas_corregidas_var2))
## [1] 144
print(length(tasa_crecimiento_var2))
## [1] 144
print(length(tasa_tendencia_var2))
## [1] 144

Gráfico variable original y tendencia variable 2, AZÚCAR (Producción de azúcar): tasa de crecimiento anual

# 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 = "#EEA2AD", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable 2 - AZÚCAR (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_var2)

Gráfico variable original y tendencia variable 3, X_AZU (Exportación de azúcar): tasa de crecimiento anual

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var3 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var3))

# Verificar longitudes
print(length(fechas_corregidas_var3))
## [1] 144
print(length(tasa_crecimiento_var3))
## [1] 144
print(length(tasa_tendencia_var3))
## [1] 144

Gráfico variable original y tendencia variable 3, X_AZU (Exportación de azúcar): tasa de crecimiento anual

# 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 = "#EEA2AD", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable 3 - X_AZU (Exportació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_var3)

La molienda de caña es el proceso fundamental y el insumo primario para la producción de azúcar. Existe una relación de causa y efecto directa, sin molienda de caña, no hay producción de azúcar ni exportación de la misma. Esto la convierte en el principal indicador adelantado de la producción.Como se observó en las gráficas, la variable de molienda se correlaciona casi perfectamente con la variable de producción de azúcar. La estacionalidad y los picos de una variable se replican fielmente en la otra. Esta relación tan estrecha y predecible la hace ideal para construir un modelo de predicción robusto.Las gráficas confirman la estrecha relación entre los procesos de la cadena de valor del azúcar: la molienda es el motor principal de la producción, y la producción es el principal impulsor de las exportaciones.

Modelo Arima y Sarima

Modelo ARIMA para la variable 1, CAN (Molienda de azúcar)

# 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

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(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2        mean
##       -0.8309  1.4498  0.6286  1916552.00
## s.e.   0.0837  0.0785  0.0713    45237.36
## 
## sigma^2 = 1.143e+11:  log likelihood = -2163.34
## AIC=4336.68   AICc=4337.09   BIC=4351.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 441.7723 333667.1 262674.3 -9.998612 22.11749 1.401977 -0.01737076

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       -8.3086e-01  8.3676e-02 -9.9295 < 2.2e-16 ***
## ma1        1.4498e+00  7.8479e-02 18.4733 < 2.2e-16 ***
## ma2        6.2862e-01  7.1328e-02  8.8130 < 2.2e-16 ***
## intercept  1.9166e+06  4.5237e+04 42.3666 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
darima_auto <- Arima(train_ts, 
                order = c(1, 0, 2))  # Especificamos directamente (p=4, d=1, q=2)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2        mean
##       -0.8309  1.4498  0.6286  1916552.00
## s.e.   0.0837  0.0785  0.0713    45237.36
## 
## sigma^2 = 1.143e+11:  log likelihood = -2163.34
## AIC=4336.68   AICc=4337.09   BIC=4351.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 441.7723 333667.1 262674.3 -9.998612 22.11749 1.401977 -0.01737076

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(1,0,2) with non-zero mean
## Q* = 175.49, 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("Toneladas")

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

Interpretación

El gráfico muestra una clara discrepancia entre las toneladas de azúcar pronosticadas (línea azul) y las observadas (línea roja).

Pronóstico (Línea Azul): El modelo pronostica una leve disminución al inicio, seguida de una estabilización en la producción. Esto sugiere un escenario de molienda constante y predecible, sin mayores contratiempos.

Observado (Línea Roja): La realidad fue diferente, pues las toneladas de azúcar molidas cayeron drásticamente, lo que indica una interrupción significativa y abrupta en las operaciones. Posteriormente, la producción se recuperó con una subida fuerte y rápida, lo que podría interpretarse como un esfuerzo para compensar la caída inicial.

En resumen, mientras el modelo preveía un proceso de molienda casi ininterrumpido, la operación real enfrentó un fuerte revés, seguido de un esfuerzo considerable por recuperar el ritmo perdido. Esto subraya la necesidad de mejorar el modelo de pronóstico para que sea más sensible a los factores que pueden causar interrupciones en la producción.

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   2203426      2096791
## 2 2024.833   1856744      1915757
## 3 2024.917   2082758      1917212

Identificación automática del modelo SARIMA para la variable 1, CAN (Molienda de azúcar)

# 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(2,0,0)(2,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2     sma1      drift
##       0.3062  0.3471  -0.2317  -0.1640  -0.6308  -198.5052
## s.e.  0.0813  0.0795   0.1466   0.1217   0.1371  1201.0995
## 
## sigma^2 = 3.591e+10:  log likelihood = -1916.58
## AIC=3847.16   AICc=3848.01   BIC=3867.81

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

# Cargar el paquete necesario
library(forecast)

darima <- Arima(train_ts, 
                order = c(2, 0, 0),  
                seasonal = list(order = c(2, 1, 1), 
                                period = 12))  # Periodicidad estacional de 12 meses
summary(darima)
## Series: train_ts 
## ARIMA(2,0,0)(2,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2     sma1
##       0.3066  0.3477  -0.2323  -0.1641  -0.6300
## s.e.  0.0813  0.0794   0.1466   0.1218   0.1372
## 
## sigma^2 = 3.565e+10:  log likelihood = -1916.6
## AIC=3845.19   AICc=3845.82   BIC=3862.88
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 6199.691 178020.1 120147.6 -4.797893 11.05979 0.6412666
##                     ACF1
## Training set -0.03951852

Validación de residuales o errores del modelo

# 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(2,0,0)(2,1,1)[12]
## Q* = 37.707, df = 19, p-value = 0.006466
## 
## Model df: 5.   Total lags used: 24

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

# 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("Toneladas")

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

Interpretación

El gráfico del análisis de la molienda de azúcar revela que, el modelo de pronóstico subestima consistentemente el volumen de producción real. A pesar de seguir la misma tendencia, el pronóstico se mantiene por debajo de lo observado, generando un riesgo significativo para la gestión operativa.

Esta subestimación puede causar escasez de recursos esenciales como repuestos, energía y personal, ya que la planificación se basa en cifras incorrectas. Para asegurar la eficiencia y rentabilidad, es crucial ajustar y calibrar el modelo de pronóstico, permitiendo una asignación de recursos más precisa y una mejor toma de decisiones.

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   2203426      2139263
## 2 2024.833   1856744      1578922
## 3 2024.917   2082758      2012364

Resultados

El primer análisis mostró una clara discrepancia entre el pronóstico y la producción real de azúcar. El modelo no anticipó una caída abrupta en la molienda, lo que llevó a una baja repentina seguida de un esfuerzo de recuperación. Esta falla subraya la importancia de un modelo de pronóstico que sea más sensible a las interrupciones inesperadas.

El segundo análisis identificó que el modelo subestimó consistentemente las toneladas de azúcar molidas. Este error puede generar problemas operativos como la escasez de recursos, un mal registro de la producción y la ineficiente planificación de ventas e inventario. Esto puede afectar directamente los compromisos con los clientes si se basan en pronósticos inexactos.

Pronóstico y evaluación del modelo

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

# 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    2096791
## 2 2024.833    1915757
## 3 2024.917    1917212
## 4 2025.000    1916003
# 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 = 1916003.29658617"

Pronóstico del modelo automático SARIMA fuera de muestra, es decir, en enero 2025

# 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    2139263
## 2 2024.833    1578922
## 3 2024.917    2012364
## 4 2025.000    1941393
# 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 = 1941393.1617685"

Interpretación

El modelo ARIMA pronostica una producción de 1,916,003 toneladas para enero de 2025. Esta predicción representa un crecimiento anual del 19.59% en comparación con enero de 2024 (1.602.184 toneladas), pero un descenso mensual del -8.01% respecto a diciembre de 2024 (2.082.758 toneladas).

El modelo SARIMA pronostica 1.941.393 toneladas para el mismo período. Esto indica un crecimiento anual del 21.17% desde enero de 2024 y una disminución mensual del -6.79% desde diciembre de 2024.

Ambos modelos anticipan correctamente una desaceleración en la producción de diciembre a enero. Sin embargo, el modelo SARIMA predice un tonelaje general más alto para enero de 2025 que el modelo ARIMA (1.941.393 frente a 1.916.003).

El pronóstico más alto del modelo SARIMA y su menor disminución mensual (-6.79% frente a -8.01%) sugieren que este modelo captura las tendencias estacionales de manera más efectiva. Dado que los modelos SARIMA están diseñados específicamente para manejar la estacionalidad, este resultado es coherente. El crecimiento anual del 21.17% pronosticado por el modelo SARIMA es un indicador muy fuerte de que el sector se está expandiendo significativamente, a pesar de las fluctuaciones típicas de mes a mes. Esto sugiere una demanda robusta y el potencial para un año muy productivo.

Recomendaciones Estratégicas

📊 Se debe invertir en modernización y expansión de las instalaciones de molienda y producción para aprovechar el crecimiento del sector. Esto aumentará la capacidad y eficiencia, asegurando que la empresa se mantenga competitiva y pueda satisfacer la creciente demanda a largo plazo.

📊 La gestión estratégica de inventarios es fundamental para optimizar costos y eficiencia. Se deben usar los pronósticos para asegurar reservas suficientes durante periodos de menor producción. Esto requiere combinar una visión de corto plazo con la de largo plazo para planificar turnos, mantenimiento y compra de insumos de manera efectiva.

📊 Para un crecimiento sostenido, es crucial diversificar los mercados de exportación y analizar los mercados emergentes. También se debe aprovechar el sólido desempeño anual como herramienta de marketing para atraer nuevos inversionistas. Además, es estratégico para una empresa del sector producir azúcar para generar etanol y energía. Esta práctica no solo diversifica el negocio al acceder al mercado de combustibles renovables, sino que también permite aprovechar subproductos como el bagazo para la cogeneración de energía, reduciendo la dependencia de los combustibles fósiles. Esto, a su vez, genera nuevos flujos de ingresos y mejora la sostenibilidad ambiental de la operación.

📊 El análisis de la estacionalidad es clave para entender las variaciones de la producción, lo que permite una mejor planificación de recursos y la anticipación de paros climáticos. Para mejorar aún más la precisión del pronóstico, es crucial ajustarlo por errores recurrentes e integrar variables externas como la lluvia y la temperatura, lo que resulta en predicciones más fiables basadas en señales que influyen directamente en la molienda.