CONTEXTO DEL SECTOR Y EMPRESA

1- Sector económico: Agroindrustrial 2- Empresa: Riopaila Castilla

Elegimos la empresa Riopaila Castilla porque es un grupo que pertenece al sector agroindustrial colombiano con más de 100 años de experiencia en transformar la tierra en alimento (azúcar, mieles, jarabes, mezclas, aceite de palma y ganadería), energía verde (cogeneración) y combustible renovable (alcohol carburante).

Riopaila Castilla opera en 36 municipios de Colombia, ubicados en los departamentos del Valle del Cauca, Cauca y Vichada. Con sus productos, llega a más de 15 millones de hogares, cogenera energía para iluminar a más de 110 mil personas y aporta a la oxigenación, en promedio, de 440 mil vehículos en el país con combustible renovable. Además, tiene presencia en 46 países con sus productos.

SELECCION Y JUSTIFICACION DE LAS 3 VARIABLES:

3- Variables relevantes: Producción de azúcar, Molienda de caña de azúcar y Exportaciones de azúcares y confitería

Se seleccionan estas tres variables porque permiten analizar al ingenio Riopaila Castilla desde sus tres dimensiones principales: producción, operación y mercado. la produccion de azucar mide el desempeño productivo del ingenio, la molienda de caña de azucar refleja eficiencia del molino y la continuidad operativa, y exportaciones de azucares y confiteria incorpora el comportamiento comercial asociado a exportaciones del sector azucarero.

Para las tres variables seleccionadas se esperaría una relación complementaria, las variables de produccion y molienda deberían comportarse de forma relacionada, porque la molienda refleja la capacidad y eficiencia operativa del ingenio, mientras que la produccion representa el resultado productivo final. Por otro lado la exportacion añadiría la dimensión comercial: captaría cómo esa producción se traduce en exportaciones del sector, por lo que su comportamiento puede acompañar a la produccion, aunque también depender de condiciones de mercado y política comercial.

Esta información puede apoyar decisiones de planificación de producción e inventarios en Riopaila Castilla, porque permite anticipar si la capacidad de molienda y la producción de azúcar estarán alineadas con la demanda comercial/exportadora. También puede servir para ajustar la operación del ingenio y decidir cuánto producir, cuánto almacenar y cuánto destinar a exportación según el comportamiento esperado de las variables.

METOLODGIA

Instalar/Cargar librerias necesarias para el análisis

#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("C:/Users/fabia/OneDrive/Desktop/MAESTRIA/ANALITICA DE NEGOCIOS/MODULO 2/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"))

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

Variable 1

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

Variable 2

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

Variable 3

# Convertir/declarar Exportaciones de azúcares y confitería  en serie de tiempo mensual
variable3_ts <- ts(data_col$X_AZU, start = c(2012, 1), frequency = 12)

EXTRACCION DE SEÑALES

Gráfico inicial de la variable 1 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("Figura 1- Variable 1: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

En la figura 1 se observa que la serie original de producción de azúcar presenta una alta variabilidad a lo largo del periodo analizado, con fluctuaciones frecuentes entre valores altos y bajos. A simple vista no se identifica una tendencia creciente o decreciente marcada, aunque se observan algunos picos y caídas pronunciadas. Asimismo, se aprecia cierta repetición en el comportamiento de la serie, lo que podría indicar la presencia de un componente estacional.

Extracción señales variable 1

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

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

# Crear fechas reales mensuales
fechas <- seq.Date(
  from = as.Date("2012-01-01"),
  by = "month",
  length.out = length(variable1_ts)
)

# Crear dataframe
stl_df_var1 <- data.frame(
  
  Time = rep(fechas, 4),
  
  Value = c(
    stl_decomp_var1$time.series[, "seasonal"],
    stl_decomp_var1$time.series[, "trend"],
    stl_decomp_var1$time.series[, "remainder"],
    as.numeric(variable1_ts)
  ),
  
  Component = rep(
    c(
      "Estacional",
      "Tendencia",
      "Residuo",
      "Serie Original"
    ),
    each = length(variable1_ts)
  )
)

# Gráfico
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 = "Figura 2 - Descomposición temporal de la variable 1",
    x = "Tiempo",
    y = "Valor"
  )

# Convertir a gráfico interactivo
ggplotly(p)

En la figura 2 se evidencia una estacionalidad definida en la producción de azúcar, caracterizada por un patrón recurrente a lo largo de los años, donde se muestran caidas para el mes de mayo y picos para mes de agosto aproximadamente. La tendencia presenta variaciones durante el periodo, con una disminución general en los últimos años y una leve recuperación al final. Por su parte, en los residuos se identifica un valor extremo alrededor del mes de mayo de 2021 que podria ser causado el paro nacional en Colombia, iniciado a finales de abril, que generó bloqueos, disrupciones logísticas y afectaciones a la actividad productiva. En conjunto, estos resultados indican que la estacionalidad constituye un componente relevante para explicar y pronosticar el comportamiento de la producción de azucar.

Extracción señales variable 2

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

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

# Crear fechas reales mensuales
fechas_var2 <- seq.Date(
  from = as.Date("2012-01-01"),
  by = "month",
  length.out = length(variable2_ts)
)

# Crear dataframe para graficar
stl_df_var2 <- data.frame(
  
  Time = rep(fechas_var2, 4),
  
  Value = c(
    stl_decomp_var2$time.series[, "seasonal"],
    stl_decomp_var2$time.series[, "trend"],
    stl_decomp_var2$time.series[, "remainder"],
    as.numeric(variable2_ts)
  ),
  
  Component = rep(
    c(
      "Estacional",
      "Tendencia",
      "Residuo",
      "Serie Original"
    ),
    each = length(variable2_ts)
  )
)

# Crear gráfico
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 = "Figura 3- Descomposición temporal de la variable 2",
    x = "Tiempo",
    y = "Valor"
  )

# Convertir a gráfico interactivo
ggplotly(p)

En la figura 3 la descomposición temporal evidencia igualmente una estacionalidad muy marcada, con caidas repetitivas para el mes de mayo y picos en el mes de agosto, que coincide con el comportamiento estacional de la variable produccion de azucar. La tendencia presenta diferentes fases de crecimiento y disminución, alcanzando niveles elevados alrededor de 2018 y mostrando posteriormente una reducción hasta 2024, seguida de una recuperación al final del periodo. En los residuos se identifica un valor atípico negativo alrededor del mes de mayo de 2021, igual que en la variable produccion de azucar.

Extracción señales variable 3

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

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

# Crear fechas reales mensuales
fechas_var3 <- seq.Date(
  from = as.Date("2012-01-01"),
  by = "month",
  length.out = length(variable3_ts)
)

# Crear dataframe para graficar
stl_df_var3 <- data.frame(
  
  Time = rep(fechas_var3, 4),
  
  Value = c(
    stl_decomp_var3$time.series[, "seasonal"],
    stl_decomp_var3$time.series[, "trend"],
    stl_decomp_var3$time.series[, "remainder"],
    as.numeric(variable3_ts)
  ),
  
  Component = rep(
    c(
      "Estacional",
      "Tendencia",
      "Residuo",
      "Serie Original"
    ),
    each = length(variable3_ts)
  )
)

# Crear gráfico
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 = "Figura 4 - Descomposición temporal de la variable 3",
    x = "Tiempo",
    y = "Valor"
  )

# Convertir a gráfico interactivo
ggplotly(p)

En la figura 4 la descomposición temporal de las exportaciones de azúcares tambien se evidencia una estacionalidad marcada, con un patrón recurrente de picos y descensos a lo largo de los años, pero en este caso las caidas se presentan para el mes de junio y los picos en el mes de octubre de cada año. La tendencia presenta diferentes fases, destacándose una disminución entre 2015 y 2019, seguida de una recuperación a partir de 2021 y niveles relativamente altos hacia el final del periodo. El componente residual presentan algunas variaciones extraordinarias, para mayo de 2021 se observa una caída, coincidente con el periodo del paro nacional en Colombia. Sin embargo, a diferencia de las variables 1 y 2, esta disminución no representa la mayor perturbación de los residuos durante el periodo analizado. Esto sugiere que, aunque el paro pudo haber afectado las exportaciones, su impacto sobre esta variable fue relativamente menor frente a otras fluctuaciones extraordinarias registradas a lo largo de la serie.

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 Variable 1

# 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 = "grey", 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("Figura 5- Variable 1:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 1") +
  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)

En la figura 5 se observa que la grafica ajustada por estacionalidad suaviza los picos y caidas marcados en el componentente estacional. Al realizar el ajuste estacional, se elimina este componente recurrente y se obtiene una serie que permite analizar con mayor claridad la tendencia y las variaciones no estacionales. Esto facilita la identificación de cambios reales en el comportamiento de las variables y mejora la interpretación de los resultados del modelo de pronóstico.

Gráfico serie original VS ajustada Variable 2

# 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("Figura 6 - Variable 2:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 2") +
  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)

En la figura 6 tambien se observa que la grafica ajustada por estacionalidad suaviza los picos y caidas marcados en el componentente estacional, lo que facilita la identificación de cambios reales en el comportamiento de las variables y mejora la interpretación de los resultados del modelo de pronóstico

Gráfico serie original VS ajustada Variable 3

# 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("Figura 7- Variable 3:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 3") +
  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)

En la figura 7 el efecto del ajuste estacional es menos evidente visualmente que en las variables anteriores. Esto puede deberser a que, aunque la serie presenta un patrón estacional recurrente, este componente tiene una menor incidencia sobre el comportamiento de la seria original. Además, la serie presenta una mayor influencia de fluctuaciones irregulares, asociadas a factores externos al patrón estacional. Por tanto, el ajuste permite depurar la variación periódica, aunque la serie ajustada mantiene una dinámica similar a la original debido al peso de los componentes de tendencia e irregularidad.

Ahora graficamos serie original vs tendencia

Primero se debe obtener la tendencia de cada variable y luego graficarla

Tendencia Variable 1

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" = "grey", "Tendencia" = "black")) +
  ggtitle("Figura 8 - Variable 1: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  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

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("Figura 9 - Variable 2: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 2") +
  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

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("Figura 10- Variable 3: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 3") +
  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)

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia:

Tasa de crecimiento de la serie de tendencia y original para la variable 1

#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] 156
print(length(tasa_crecimiento_var1))
## [1] 156
print(length(tasa_tendencia_var1))
## [1] 156

*Gráfico variable original y tendencia variable 1: tasa de crecimiento anual**

# 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("Figura 11- Variable1: Tasa de crecimiento anual % de la serie Original
          y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  coord_cartesian(ylim = c(-20, 25)) +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

En la figura 11 se evidencia que la tasa de crecimiento de la serie original presenta una volatilidad elevada, con cambios frecuentes entre periodos de crecimiento y contracción. En contraste, la tasa de crecimiento de la tendencia permite identificar de manera más clara los ciclos de expansión y desaceleración de la variable, al reducir el efecto de las fluctuaciones de corto plazo. Por ello, la tendencia resulta más útil para evaluar la evolución estructural de la variable y como referencia para el análisis y pronóstico de su comportamiento futuro.

Al final del periodo analizado, la tasa de crecimiento de la serie original continúa mostrando una volatilidad elevada ; sin embargo, la tasa de crecimiento de la tendencia permanece en terreno positivo. Esto indica que, aunque la variable presenta fluctuaciones importantes en el corto plazo, su comportamiento estructural más reciente continúa siendo creciente, aunque evidencia una desaceleración para los 3 primeros trimestres del año 2025, los últimos periodos sugieren una posible recuperación en su ritmo de crecimiento.

# Encontrar el crecimiento máximo
indice_max <- which.max(tasa_crecimiento_var1)

# Mostrar los valores que lo generan
print(indice_max)
## [1] 113
print(variable1_ts[indice_max])
## [1] 11041.41
print(variable1_ts[indice_max + 12])
## [1] 119768.7
print(tasa_crecimiento_var1[indice_max])
## [1] 984.7229
print(fechas_corregidas_var1[indice_max])
## [1] "2022-05-01"
# Revisar los 12 meses alrededor del periodo del pico
variable1_ts[100:125]
##  [1] 131681.81 114502.47 174447.47 203347.14 228369.29 228980.14 234065.74
##  [8] 146351.23 193674.35 177254.28 193573.94 145987.21 158408.21  11041.41
## [15] 156211.52 220751.10 233249.23 226637.25 211407.85 162581.72 202837.56
## [22] 197814.77 189215.64 180967.04 139870.52 119768.71
# Crear tabla con fechas y valores
datos_revision <- data.frame(
  fecha = seq(
    from = as.Date("2012-01-01"),
    by = "month",
    length.out = length(variable1_ts)
  ),
  valor = as.numeric(variable1_ts)
)

# Ver 2021 y 2022
datos_revision[
  datos_revision$fecha >= as.Date("2021-01-01") &
  datos_revision$fecha <= as.Date("2022-12-01"),
]
fecha valor
109 2021-01-01 177254.28
110 2021-02-01 193573.94
111 2021-03-01 145987.21
112 2021-04-01 158408.21
113 2021-05-01 11041.41
114 2021-06-01 156211.52
115 2021-07-01 220751.10
116 2021-08-01 233249.23
117 2021-09-01 226637.25
118 2021-10-01 211407.85
119 2021-11-01 162581.72
120 2021-12-01 202837.56
121 2022-01-01 197814.77
122 2022-02-01 189215.64
123 2022-03-01 180967.04
124 2022-04-01 139870.52
125 2022-05-01 119768.71
126 2022-06-01 118137.14
127 2022-07-01 189326.35
128 2022-08-01 217934.39
129 2022-09-01 225602.05
130 2022-10-01 201336.46
131 2022-11-01 136489.15
132 2022-12-01 178279.33

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2

#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] 156
print(length(tasa_crecimiento_var2))
## [1] 156
print(length(tasa_tendencia_var2))
## [1] 156
# 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("Figura 12- Variable2: Tasa de crecimiento anual % de la serie Original
          y la Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  coord_cartesian(ylim = c(-20, 25)) +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

En la figura 12 se evidencia que la variable 2 presenta un comportamiento similar al observado en la variable 1, tanto en la evolución de la tasa de crecimiento como en el comportamiento de su tendencia. En ambas series se observa un crecimiento importante al inicio del periodo, seguido de etapas de desaceleración y contracción. Hacia el final del periodo, ambas variables muestran una recuperación de su tasa de crecimiento, aunque la variable 2 alcanza niveles de crecimiento de tendencia más elevados durante 2025.

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 3

#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] 156
print(length(tasa_crecimiento_var3))
## [1] 156
print(length(tasa_tendencia_var3))
## [1] 156
# Gráfico de la tasa de crecimiento anual variable 3
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("Figura 13- Variable 3: Tasa de crecimiento anual % de la serie Original
          y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  coord_cartesian(ylim = c(-30, 40)) +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

En la figura 13 se observa que la variable 3 presenta una dinámica más fluctuante en su tasa de crecimiento, con cambios frecuentes entre periodos de expansión y contracción. Hacia el final del periodo, la variable 3 muestra un comportamiento diferente al de las variables 1 y 2. Después del crecimiento observado hacia finales del año 2024, la tasa de crecimiento de la tendencia comienza a desacelerarse y pasa nuevamente a valores negativos hacia mediados del años 2025. Al cierre de la serie, la tendencia se ubica en terreno negativo, lo que indica que las exportaciones presentan una pérdida de dinamismo y una contracción reciente.

ANALISIS INTEGRADO DE LAS VARIABLES

El análisis conjunto de las tres variables evidencia un escenario mixto para el sector y la empresa. Si bien se observan periodos de recuperación y crecimiento en algunas de las variables, el comportamiento no es homogéneo para las 3 variables.

Históricamente, las tres series presentan comportamientos cíclicos y fluctuaciones importantes. Un punto particularmente relevante corresponde a mayo de 2021, periodo en el que las variables 1 y 2 registran una caída bastante marcada, asociada al impacto del paro nacional. En el caso de la variable 3, correspondiente a las exportaciones de azúcares, también se observa una disminución en este periodo; sin embargo, la caída no constituye el punto más extremo de la serie, debido a que esta variable presenta una mayor volatilidad durante otros momentos del periodo analizado.

# Gráfico conjunto de las tasas de crecimiento anual
grafico_integracion_crecimiento <- ggplot() +
  
  geom_line(aes(x = fechas_corregidas_var1,
                y = tasa_tendencia_var1,
                color = "Variable 1"),
            linewidth = 0.9) +
  
  geom_line(aes(x = fechas_corregidas_var2,
                y = tasa_tendencia_var2,
                color = "Variable 2"),
            linewidth = 0.9) +
  
  geom_line(aes(x = fechas_corregidas_var3,
                y = tasa_tendencia_var3,
                color = "Variable 3"),
            linewidth = 0.9) +
  
  geom_hline(yintercept = 0,
             linetype = "dashed",
             color = "black",
             linewidth = 0.5) +
  
  ggtitle("Figura 14 - Tasas de crecimiento anual de las tres variables") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  
  scale_color_manual(
    values = c(
      "Variable 1" = "#1F4E79",
      "Variable 2" = "#5B9BD5",
      "Variable 3" = "#70AD47"
    )
  ) +
  
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 13),
    legend.title = element_blank(),
    legend.position = "bottom"
  )

# Convertir a gráfico interactivo
ggplotly(grafico_integracion_crecimiento)

En la figura 14 se observa que las señales presentan un comportamiento mixto. Las variables 1 y 2 muestran una dinámica similar y evidencian una recuperación hacia el final del periodo, mientras que la variable 3 presenta una mayor volatilidad y un comportamiento menos sincronizado. Al cierre del periodo, las variables 1 y 2 mantienen tasas de crecimiento positivas, en contraste con la variable 3, que vuelve a registrar una tasa negativa. Esto sugiere que, aunque existen señales de recuperación en parte del sector, persisten factores de incertidumbre en el desempeño de sus componentes, por lo que el escenario actual debe interpretarse con cautela.

Modelo ARIMA

División en conjunto de entrenamiento y prueba para la variable 1 que es la elegida para pronosticar

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

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(2025, 9))  # Entrenamiento hasta septiembre 2025
test_ts <- window(variable1_ts, start = c(2025, 10))  # Prueba inicia desde oct2025

Modelo ARIMA automático normal (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.0681  -0.4894  -0.4726  179805.398
## s.e.  0.1878   0.0932   0.2115    3444.423
## 
## sigma^2 = 1.269e+09:  log likelihood = -1961.69
## AIC=3933.38   AICc=3933.76   BIC=3948.91
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 37.14014 35190.76 28085.81 -12.88172 26.65939 1.371507 -0.01162577

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.0681e+00  1.8783e-01  5.6864 1.298e-08 ***
## ar2       -4.8943e-01  9.3167e-02 -5.2532 1.495e-07 ***
## ma1       -4.7257e-01  2.1155e-01 -2.2339   0.02549 *  
## intercept  1.7981e+05  3.4444e+03 52.2019 < 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.0681  -0.4894  -0.4726  179805.398
## s.e.  0.1878   0.0932   0.2115    3444.423
## 
## sigma^2 = 1.269e+09:  log likelihood = -1961.69
## AIC=3933.38   AICc=3933.76   BIC=3948.91
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 37.14014 35190.76 28085.81 -12.88172 26.65939 1.371507 -0.01162577

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* = 196.74, 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 automático (2,0,1) evidencia diferencias entre ambas trayectorias. Aunque el modelo logra representar parcialmente el nivel de la variable, no captura adecuadamente el punto de quiebre observado al final del período de evaluación. Mientras la serie observada presenta una recuperación después de una caída, el pronóstico mantiene una trayectoria descendente. Esto indica que el modelo presenta limitaciones para anticipar puntos de quiebre o cambios estructurales de corto plazo.

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),
  Observado = as.numeric(test_ts),
  Pronosticado = as.numeric(arima_forecast_auto$mean)
)

# Convertir el tiempo a formato mes-año
forecast_table_auto$Tiempo <- format(
  as.Date(as.yearmon(forecast_table_auto$Tiempo)),
  "%B %Y"
)

# Mostrar la tabla
print(forecast_table_auto)
##           Tiempo Observado Pronosticado
## 1   octubre 2025  185732.9     194842.5
## 2 noviembre 2025  161012.2     177370.0
## 3 diciembre 2025  180947.8     169844.7

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

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

# Número de períodos del conjunto de prueba
n_test <- length(test_ts)

# Pronosticar todo el período de prueba + enero, febrero y marzo de 2026
arima_forecast_2026 <- forecast(
  auto_arima_model_no_seasonal,
  h = n_test + 3
)

# Extraer únicamente los 3 meses adicionales
pronostico_2026 <- tail(
  as.numeric(arima_forecast_2026$mean),
  3
)

# Crear tabla con los tres meses
pronostico_2026_table <- data.frame(
  Mes = c("Enero 2026", "Febrero 2026", "Marzo 2026"),
  Pronostico = round(pronostico_2026, 2)
)

# Mostrar los tres pronósticos
print(pronostico_2026_table)
##            Mes Pronostico
## 1   Enero 2026   170358.7
## 2 Febrero 2026   174590.8
## 3   Marzo 2026   178859.4
library(zoo)

# Fechas de la serie de la variable 1
fechas_var1 <- as.Date(as.yearmon(time(variable1_ts)))

# Extraer los valores reales de enero, febrero y marzo de 2025
valores_2025 <- as.numeric(
  variable1_ts[
    fechas_var1 %in% as.Date(c(
      "2025-01-01",
      "2025-02-01",
      "2025-03-01"
    ))
  ]
)

# Valores pronosticados para enero, febrero y marzo de 2026
valores_2026 <- pronostico_2026

# Calcular crecimiento anual
tasa_crecimiento_2026 <- 
  ((valores_2026 / valores_2025) - 1) * 100

# Crear tabla
tabla_crecimiento_2026 <- data.frame(
  Mes = c("Enero 2026", "Febrero 2026", "Marzo 2026"),
  Valor_2025 = round(valores_2025, 2),
  Pronostico_2026 = round(valores_2026, 2),
  Tasa_crecimiento_anual = round(tasa_crecimiento_2026, 2)
)

print(tabla_crecimiento_2026)
##            Mes Valor_2025 Pronostico_2026 Tasa_crecimiento_anual
## 1   Enero 2026   170559.3        170358.7                  -0.12
## 2 Febrero 2026   153275.8        174590.8                  13.91
## 3   Marzo 2026   144025.2        178859.4                  24.19

Modelo SARIMA automático

Identificación dautomá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,0)(1,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2     ar3     sar1     sma1      drift
##       0.3328  0.2633  0.0547  -0.1548  -0.8380  -128.1560
## s.e.  0.0845  0.0826  0.0857   0.1021   0.0996    82.9244
## 
## sigma^2 = 335634360:  log likelihood = -1724.85
## AIC=3463.71   AICc=3464.48   BIC=3484.92

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(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(3, 0, 0),  # (p,d,q) -> (3,0,0)
                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,0)(1,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     ar3     sar1     sma1
##       0.3384  0.2771  0.0770  -0.1685  -0.8044
## s.e.  0.0850  0.0822  0.0836   0.1017   0.0895
## 
## sigma^2 = 341478895:  log likelihood = -1725.81
## AIC=3463.62   AICc=3464.2   BIC=3481.8
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -803.5598 17501.32 12264.55 -5.948787 12.18748 0.5989111
##                    ACF1
## Training set 0.01137898

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

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(
  auto_arima_model,
  h = length(test_ts)
)

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

# Gráfico
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico"), linewidth = 1) +
  geom_line(aes(y = Observado, color = "Observado"), linewidth = 1) +
  ggtitle("Pronóstico vs Observado - Variable 1") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(p4)

Interpretación modelo SARIMA (3,0,0)(1,1,1): Este modelo parece pronosticar mejor dentro de prueba. Hay una sub estimacion al final del periodo de prueba, pero se capturan muy bien los puntos de quiebre. Es un modelo tentativo adecuado para pronpostico fuera de muestra o a futuro.

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

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

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

# Crear fechas del periodo pronosticado
fechas <- as.Date(as.yearmon(time(arima_forecast$mean)))

# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
  Mes = format(fechas, "%B %Y"),
  Observado = as.numeric(test_ts),
  Pronosticado = as.numeric(arima_forecast$mean)
)

# Mostrar la tabla
print(forecast_table)
##              Mes Observado Pronosticado
## 1   octubre 2025  185732.9     193401.1
## 2 noviembre 2025  161012.2     132674.1
## 3 diciembre 2025  180947.8     172438.6

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

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

# Hacer el pronóstico para el conjunto de prueba + 3 meses de 2026
next_forecast <- forecast(
  auto_arima_model,
  h = length(test_ts) + 3
)

# Obtener las fechas
fechas <- as.Date(as.yearmon(time(next_forecast$mean)))

# Extraer únicamente enero, febrero y marzo de 2026
pronostico_2026 <- data.frame(
  Mes = c("Enero 2026", "Febrero 2026", "Marzo 2026"),
  Pronostico = round(
    as.numeric(next_forecast$mean[
      fechas %in% as.Date(c(
        "2026-01-01",
        "2026-02-01",
        "2026-03-01"
      ))
    ]), 
    2
  )
)

# Mostrar únicamente los tres meses
print(pronostico_2026)
##            Mes Pronostico
## 1   Enero 2026   154286.4
## 2 Febrero 2026   174984.9
## 3   Marzo 2026   172137.4
# Valores reales de enero, febrero y marzo de 2025
fechas_var1 <- as.Date(as.yearmon(time(variable1_ts)))

valores_2025 <- as.numeric(
  variable1_ts[
    fechas_var1 %in% as.Date(c(
      "2025-01-01",
      "2025-02-01",
      "2025-03-01"
    ))
  ]
)

# Valores pronosticados de enero, febrero y marzo de 2026
valores_2026 <- pronostico_2026$Pronostico

# Calcular tasa de crecimiento anual
tasa_crecimiento <- ((valores_2026 / valores_2025) - 1) * 100

# Crear tabla final
tabla_crecimiento_2026 <- data.frame(
  Mes = c("Enero 2026", "Febrero 2026", "Marzo 2026"),
  Valor_2025 = round(valores_2025, 2),
  Pronostico_2026 = round(valores_2026, 2),
  Tasa_crecimiento_anual = round(tasa_crecimiento, 2)
)

# Mostrar resultado
print(tabla_crecimiento_2026)
##            Mes Valor_2025 Pronostico_2026 Tasa_crecimiento_anual
## 1   Enero 2026   170559.3        154286.4                  -9.54
## 2 Febrero 2026   153275.8        174984.9                  14.16
## 3   Marzo 2026   144025.2        172137.4                  19.52

CONCLUSIONES

RECOMENDACIONES