Resumen ejecutivo

  • Señal. La tendencia-ciclo de los despachos de cemento al Valle del Cauca crece −3,7 % anual en diciembre de 2025, mientras el mercado nacional crece +9,2 %: el problema es regional, no del país. El área licenciada en Cali cae −36,6 % en 2025 y su tendencia va −77,7 % anual, el peor registro de la serie.
  • Pronóstico. El modelo seleccionado, un SARIMA(0,1,4)(0,0,1)[12], proyecta para enero–marzo de 2026 cerca de 260.000 t (+3,7 % anual).
  • Riesgo. El dato real publicado después por el DANE fue −2,6 % en el trimestre: el modelo sobreestimó la demanda. Extrapolar la inercia de la serie sin incorporar la señal de las licencias lleva a sobreproducir.
  • Decisión. Planear producción e inventarios en el borde inferior del intervalo de predicción, no en la media; reasignar esfuerzo comercial del canal granel (obra formal) al canal empacado (ferretería y remodelación).
  • Indicador a monitorear. Área licenciada en Cali (DANE-ELIC): si su variación anual sigue por debajo de −25 % durante tres meses más, el deterioro de los despachos se prolonga hacia 2027.

1. Selección: sector, empresa y variables

Sector. Construcción, específicamente la cadena del cemento. Es un sector cíclico, intensivo en capital, con demanda concentrada geográficamente y con un encadenamiento conocido —licencia → preventa → inicio de obra → compra de cemento— que permite construir un indicador anticipado.

Empresa. Cementos Argos S.A., operación Suroccidente, con planta en Yumbo (Valle del Cauca). Los despachos con destino al Valle son la mejor aproximación pública a la demanda que atiende esa planta. La capacidad instalada y los volúmenes propios de la planta de Yumbo no se verificaron con fuentes públicas de la compañía; por esa razón no se citan cifras propias de la planta.

Variables seleccionadas (Base Caso 1, hoja Caso2; mensuales, enero 2012 – diciembre 2025):

Código Descripción Unidad Fuente Rol en el análisis
CEM_V Despachos de cemento gris con destino al Valle del Cauca toneladas DANE – ECG Variable a pronosticar (variable 1)
LICC_CALI Área aprobada en licencias de construcción, Cali m² DANE – ELIC Indicador líder (variable 2)
DECEM Despachos nacionales de cemento gris toneladas DANE – ECG Referencia nacional (variable 3)
LICC Área licenciada nacional m² DANE – ELIC Auxiliar, solo para imputar un dato faltante

Justificación.

  1. CEM_V es la demanda del mercado que abastece la planta de Yumbo: pronosticarla equivale a pronosticar las toneladas que la operación debe producir, despachar y financiar.
  2. LICC_CALI mide el área que las curadurías aprobaron hoy y que se convertirá en obra —y por tanto en compra de cemento— varios meses después. Es la variable que puede dar aviso temprano.
  3. DECEM permite separar lo regional de lo nacional: si el Valle cae mientras Colombia crece, la causa es local y la respuesta comercial también debe serlo.

Instalar/Cargar librerías necesarias para el análisis

# Cargar librerías necesarias
library(readxl)    # Para leer archivos Excel
library(tseries)   # Para pruebas de estacionariedad (ADF, KPSS)
library(forecast)  # Para modelado ARIMA, validación cruzada y pronósticos
library(lmtest)    # Para la significancia de los coeficientes y la prueba de Granger
library(sandwich)  # Para errores estándar robustos (Newey-West)
library(ggplot2)   # Para visualización de datos
library(plotly)    # Para gráficos interactivos

Cargar base de datos

library(readxl)

# La base debe estar en la misma carpeta que este archivo (Session > Set Working Directory)
data_col <- read_excel("Base Caso1.xlsx", sheet = "Caso2")

# Nos quedamos con la fecha y las cuatro variables de interés
data_col <- data_col[, c("FECHA", "CEM_V", "LICC_CALI", "DECEM", "LICC")]
data_col$FECHA <- as.Date(data_col$FECHA)

cat("Periodo:", format(min(data_col$FECHA), "%Y-%m"), "a", format(max(data_col$FECHA), "%Y-%m"),
    "| observaciones:", nrow(data_col), "\n")
## Periodo: 2012-01 a 2025-12 | observaciones: 168
colSums(is.na(data_col))
##     FECHA     CEM_V LICC_CALI     DECEM      LICC 
##         0         0         1         0         0

Tratamiento del único dato faltante

LICC_CALI no tiene dato en abril de 2020, el mes del confinamiento nacional. Borrarlo o interpolarlo suavizaría artificialmente el choque. Se imputa con la participación mediana de Cali en el área licenciada nacional de los 12 meses previos, aplicada al dato nacional de ese mismo mes: así se conserva la caída sistémica de abril de 2020.

# Imputación del faltante de abril 2020 (participación mediana de Cali en el total nacional)
i   <- which(is.na(data_col$LICC_CALI))   # posición del faltante
idx <- (i - 12):(i - 1)                   # los 12 meses previos

participacion <- median(data_col$LICC_CALI[idx] / data_col$LICC[idx])
data_col$LICC_CALI[i] <- data_col$LICC[i] * participacion

cat("Mes imputado:", format(data_col$FECHA[i], "%Y-%m"),
    "| participación mediana de Cali:", round(100 * participacion, 2), "%",
    "| valor imputado:", round(data_col$LICC_CALI[i]), "m²\n")
## Mes imputado: 2020-04 | participación mediana de Cali: 3.9 % | valor imputado: 14539 m²

Declaración de las variables como series de tiempo

Variable 1

# Convertir/declarar variable 1 = CEM_V (despachos de cemento al Valle) en serie de tiempo mensual
variable1_ts <- ts(data_col$CEM_V, start = c(2012, 1), frequency = 12)
window(variable1_ts, start = c(2025, 1))   # últimos 12 meses
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2025 77510.19 83074.90 90124.49 84565.81 94011.39 80795.56 99417.04 90647.74
##           Sep      Oct      Nov      Dec
## 2025 95562.29 94707.63 88954.69 80316.50

Variable 2

# Convertir/declarar variable 2 = LICC_CALI (licencias de construcción en Cali) en serie de tiempo mensual
variable2_ts <- ts(data_col$LICC_CALI, start = c(2012, 1), frequency = 12)
window(variable2_ts, start = c(2025, 1))
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2025 158992  75754  35326  35455  42353  68818  87840  68070  24777  14257
##         Nov    Dec
## 2025  16177  19210

Variable 3

# Convertir/declarar variable 3 = DECEM (despachos nacionales de cemento) en serie de tiempo mensual
variable3_ts <- ts(data_col$DECEM, start = c(2012, 1), frequency = 12)
window(variable3_ts, start = c(2025, 1))
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2025  889613.5  981553.2 1084450.0 1012146.3 1099239.1  974792.5 1203723.8
##            Aug       Sep       Oct       Nov       Dec
## 2025 1098571.4 1181978.7 1201023.9 1112120.4 1071996.2

2. Estadística descriptiva

# Tabla descriptiva de las tres variables
descriptivos <- function(x, nombre, unidad) {
  anual <- tapply(x, format(data_col$FECHA, "%Y"), sum)
  data.frame(
    Variable   = nombre,
    Unidad     = unidad,
    Media      = round(mean(x)),
    Desv_est   = round(sd(x)),
    CV_pct     = round(100 * sd(x) / mean(x), 1),
    Minimo     = round(min(x)),
    Mes_min    = format(data_col$FECHA[which.min(x)], "%Y-%m"),
    Maximo     = round(max(x)),
    Mes_max    = format(data_col$FECHA[which.max(x)], "%Y-%m"),
    Total_2024 = round(anual[["2024"]]),
    Total_2025 = round(anual[["2025"]]),
    Var_2025_pct = round(100 * (anual[["2025"]] / anual[["2024"]] - 1), 1),
    TCAC_pct   = round(100 * ((anual[["2025"]] / anual[["2012"]])^(1 / 13) - 1), 2)
  )
}

tabla_desc <- rbind(
  descriptivos(data_col$CEM_V,     "Cemento Valle",    "toneladas"),
  descriptivos(data_col$LICC_CALI, "Licencias Cali",   "m² aprobados"),
  descriptivos(data_col$DECEM,     "Cemento Colombia", "toneladas")
)
tabla_desc
Variable Unidad Media Desv_est CV_pct Minimo Mes_min Maximo Mes_max Total_2024 Total_2025 Var_2025_pct TCAC_pct
Cemento Valle toneladas 84963 13814 16.3 11485 2021-05 119923 2022-09 1076366 1059688 -1.5 1.89
Licencias Cali m² aprobados 83017 54708 65.9 6524 2023-02 326535 2024-12 1020139 647029 -36.6 -3.56
Cemento Colombia toneladas 1015369 116692 11.5 242414 2020-04 1257125 2022-03 12247170 12911209 5.4 1.61 
# Participación del Valle en los despachos nacionales
participacion_valle <- data.frame(
  Anio = as.numeric(names(tapply(data_col$CEM_V, format(data_col$FECHA, "%Y"), sum))),
  Participacion_pct = round(100 * tapply(data_col$CEM_V, format(data_col$FECHA, "%Y"), sum) /
                                  tapply(data_col$DECEM, format(data_col$FECHA, "%Y"), sum), 2)
)
tail(participacion_valle, 6)
Anio Participacion_pct
2020 2020 8.90
2021 2021 8.17
2022 2022 9.24
2023 2023 8.53
2024 2024 8.79
2025 2025 8.21

Lectura. El Valle pesa alrededor del 8–9 % de los despachos nacionales, pero esa participación pasó de 8,79 % en 2024 a 8,21 % en 2025: el mercado regional no solo cae, sino que pierde peso relativo dentro de Colombia. La variable a pronosticar tiene un coeficiente de variación de 16,3 %, y su mínimo histórico (11.485 t, mayo de 2021) corresponde al Paro Nacional, no a un fenómeno de demanda.

Gráfico inicial de la variable 1 en niveles - Original

library(ggplot2)
library(plotly)

grafico_serie <- ggplot(data_col, aes(x = FECHA, y = CEM_V)) +
  geom_line(color = "grey", linewidth = 0.4) +
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Despachos de cemento al Valle del Cauca (serie original)") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  theme_minimal()

ggplotly(grafico_serie)

Gráfico inicial de la variable 2 en niveles - Original

grafico_serie2 <- ggplot(data_col, aes(x = FECHA, y = LICC_CALI)) +
  geom_line(color = "grey", linewidth = 0.4) +
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 2: Área aprobada en licencias de construcción, Cali (serie original)") +
  xlab("Tiempo") +
  ylab("m² aprobados") +
  theme_minimal()

ggplotly(grafico_serie2)

Gráfico inicial de la variable 3 en niveles - Original

grafico_serie3 <- ggplot(data_col, aes(x = FECHA, y = DECEM)) +
  geom_line(color = "grey", linewidth = 0.4) +
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 3: Despachos nacionales de cemento gris (serie original)") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  theme_minimal()

ggplotly(grafico_serie3)

Extracción de señales

Muchas series de tiempo son una combinación de varias influencias. Separar la tendencia, la estacionalidad y el componente irregular permite entender mejor qué está impulsando los cambios en la serie.

En este caso concreto, la pregunta de negocio es: ¿los despachos al Valle caen porque el mercado se está contrayendo (tendencia), porque enero y febrero son meses flojos (estacionalidad) o porque ocurrió un hecho puntual (irregular)? Cada respuesta implica una decisión distinta: la primera obliga a revisar capacidad, la segunda a programar mantenimientos, y la tercera a no hacer nada estructural.

Además, los modelos de pronóstico funcionan mejor cuando las señales subyacentes están bien definidas, y detectar cambios inesperados es más fácil cuando se eliminan los componentes predecibles.

Dos decisiones metodológicas antes de descomponer:

  1. Se descompone el logaritmo de cada serie, no el nivel. Esto equivale a una descomposición multiplicativa (Serie = Tendencia × Estacional × Irregular), que es la adecuada cuando la amplitud de las oscilaciones crece con el nivel de la serie —como ocurre en series de volumen— y permite expresar la estacionalidad en porcentaje, que es como se interpreta en la práctica. En las licencias de Cali, además, evita que la tendencia estimada tome valores negativos, algo imposible para un área licenciada.
  2. Se usa robust = TRUE en stl(). Sin esa opción, las caídas de abril de 2020 (−75 %) y mayo de 2021 (−87 %) se reparten entre el factor estacional de esos meses y contaminan el patrón de todos los años. Con la versión robusta, la descomposición asigna esos meses al componente irregular, que es donde corresponde.

Extracción de señales variable 1

# Descomposición STL sobre el logaritmo (descomposición multiplicativa) y robusta a los choques
stl_decomp_var1 <- stl(log(variable1_ts), s.window = "periodic", robust = TRUE)

# Se recuperan los componentes en unidades interpretables
tendencia_var1  <- exp(stl_decomp_var1$time.series[, "trend"])                     # toneladas
estacional_var1 <- 100 * (exp(stl_decomp_var1$time.series[, "seasonal"]) - 1)      # % frente a la tendencia
irregular_var1  <- 100 * (exp(stl_decomp_var1$time.series[, "remainder"]) - 1)     # % no explicado
variable1_sa    <- exp(log(variable1_ts) - stl_decomp_var1$time.series[, "seasonal"])  # serie ajustada

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var1 <- data.frame(
  Tiempo = rep(time(variable1_ts), 4),
  Valor  = c(as.numeric(estacional_var1), as.numeric(tendencia_var1),
             as.numeric(irregular_var1), as.numeric(variable1_ts)),
  Componente = rep(c("1. Estacional (%)", "2. Tendencia (t)", "3. Irregular (%)", "4. Serie original (t)"),
                   each = length(variable1_ts))
)

p <- ggplot(stl_df_var1, aes(x = Tiempo, y = Valor, color = Componente)) +
  geom_line() +
  facet_wrap(~ Componente, scales = "free_y", ncol = 1) +
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 1 (cemento Valle)",
       x = "Tiempo", y = "Valor")

ggplotly(p, height = 700)

Extracción de señales variable 2

stl_decomp_var2 <- stl(log(variable2_ts), s.window = "periodic", robust = TRUE)

tendencia_var2  <- exp(stl_decomp_var2$time.series[, "trend"])
estacional_var2 <- 100 * (exp(stl_decomp_var2$time.series[, "seasonal"]) - 1)
irregular_var2  <- 100 * (exp(stl_decomp_var2$time.series[, "remainder"]) - 1)
variable2_sa    <- exp(log(variable2_ts) - stl_decomp_var2$time.series[, "seasonal"])

stl_df_var2 <- data.frame(
  Tiempo = rep(time(variable2_ts), 4),
  Valor  = c(as.numeric(estacional_var2), as.numeric(tendencia_var2),
             as.numeric(irregular_var2), as.numeric(variable2_ts)),
  Componente = rep(c("1. Estacional (%)", "2. Tendencia (m²)", "3. Irregular (%)", "4. Serie original (m²)"),
                   each = length(variable2_ts))
)

p <- ggplot(stl_df_var2, aes(x = Tiempo, y = Valor, color = Componente)) +
  geom_line() +
  facet_wrap(~ Componente, scales = "free_y", ncol = 1) +
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 2 (licencias Cali)",
       x = "Tiempo", y = "Valor")

ggplotly(p, height = 700)

Extracción de señales variable 3

stl_decomp_var3 <- stl(log(variable3_ts), s.window = "periodic", robust = TRUE)

tendencia_var3  <- exp(stl_decomp_var3$time.series[, "trend"])
estacional_var3 <- 100 * (exp(stl_decomp_var3$time.series[, "seasonal"]) - 1)
irregular_var3  <- 100 * (exp(stl_decomp_var3$time.series[, "remainder"]) - 1)
variable3_sa    <- exp(log(variable3_ts) - stl_decomp_var3$time.series[, "seasonal"])

stl_df_var3 <- data.frame(
  Tiempo = rep(time(variable3_ts), 4),
  Valor  = c(as.numeric(estacional_var3), as.numeric(tendencia_var3),
             as.numeric(irregular_var3), as.numeric(variable3_ts)),
  Componente = rep(c("1. Estacional (%)", "2. Tendencia (t)", "3. Irregular (%)", "4. Serie original (t)"),
                   each = length(variable3_ts))
)

p <- ggplot(stl_df_var3, aes(x = Tiempo, y = Valor, color = Componente)) +
  geom_line() +
  facet_wrap(~ Componente, scales = "free_y", ncol = 1) +
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 3 (cemento Colombia)",
       x = "Tiempo", y = "Valor")

ggplotly(p, height = 700)

Serie original vs. serie ajustada por estacionalidad

Después de la descomposición temporal, se grafica la serie ajustada por estacionalidad junto con la original. La serie ajustada es la que debe mirar la gerencia mes a mes: responde a la pregunta “¿este mes fue bueno descontando que siempre lo es?”.

Gráfico serie original VS ajustada Variable 1

fechas <- data_col$FECHA   # vector de fechas común a las tres series

grafico_ajustada_var1 <- ggplot() +
  geom_line(aes(x = fechas, y = as.numeric(variable1_ts), color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas, y = as.numeric(variable1_sa),  color = "Serie ajustada"), linewidth = 0.6) +
  scale_color_manual(values = c("Serie original" = "grey", "Serie ajustada" = "black")) +
  ggtitle("Variable 1: Serie original vs. serie ajustada por estacionalidad") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(grafico_ajustada_var1)

Gráfico serie original VS ajustada Variable 2

grafico_ajustada_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = as.numeric(variable2_ts), color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas, y = as.numeric(variable2_sa),  color = "Serie ajustada"), linewidth = 0.6) +
  scale_color_manual(values = c("Serie original" = "grey", "Serie ajustada" = "black")) +
  ggtitle("Variable 2: Serie original vs. serie ajustada por estacionalidad") +
  xlab("Tiempo") + ylab("m² aprobados") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(grafico_ajustada_var2)

Gráfico serie original VS ajustada Variable 3

grafico_ajustada_var3 <- ggplot() +
  geom_line(aes(x = fechas, y = as.numeric(variable3_ts), color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas, y = as.numeric(variable3_sa),  color = "Serie ajustada"), linewidth = 0.6) +
  scale_color_manual(values = c("Serie original" = "grey", "Serie ajustada" = "black")) +
  ggtitle("Variable 3: Serie original vs. serie ajustada por estacionalidad") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(grafico_ajustada_var3)

Serie original vs. tendencia

La extracción de la tendencia permite centrarse en los cambios estructurales de la serie y anticipar escenarios futuros.

Tendencia Variable 1

grafico_tendencia_var1 <- ggplot() +
  geom_line(aes(x = fechas, y = as.numeric(variable1_ts),  color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas, y = as.numeric(tendencia_var1), color = "Tendencia"), linewidth = 0.9) +
  scale_color_manual(values = c("Serie original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 1: Serie original vs. tendencia") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(grafico_tendencia_var1)

Tendencia Variable 2

grafico_tendencia_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = as.numeric(variable2_ts),  color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas, y = as.numeric(tendencia_var2), color = "Tendencia"), linewidth = 0.9) +
  scale_color_manual(values = c("Serie original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 2: Serie original vs. tendencia") +
  xlab("Tiempo") + ylab("m² aprobados") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(grafico_tendencia_var2)

Tendencia Variable 3

grafico_tendencia_var3 <- ggplot() +
  geom_line(aes(x = fechas, y = as.numeric(variable3_ts),  color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas, y = as.numeric(tendencia_var3), color = "Tendencia"), linewidth = 0.9) +
  scale_color_manual(values = c("Serie original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 3: Serie original vs. tendencia") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(grafico_tendencia_var3)

Tasa de crecimiento anual de la serie original vs. la tendencia

La tendencia se expresa en tasas de crecimiento anual: el nivel en toneladas no dice si el mercado mejora o empeora, y la variación interanual sí. La tasa calculada sobre la serie original es ruidosa; la calculada sobre la tendencia muestra la dirección real del mercado.

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 <- (as.numeric(variable1_ts)[13:length(variable1_ts)] /
                            as.numeric(variable1_ts)[1:(length(variable1_ts) - 12)] - 1) * 100
tasa_tendencia_var1   <- (as.numeric(tendencia_var1)[13:length(tendencia_var1)] /
                            as.numeric(tendencia_var1)[1:(length(tendencia_var1) - 12)] - 1) * 100

# Vector de fechas corregido: la primera tasa anual corresponde a enero de 2013
fechas_corregidas <- seq(from = as.Date("2013-01-01"), by = "month",
                         length.out = length(tasa_crecimiento_var1))

cat("Longitudes:", length(fechas_corregidas), length(tasa_crecimiento_var1), length(tasa_tendencia_var1), "\n")
## Longitudes: 156 156 156

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

grafico_crecimiento_var1 <- ggplot() +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_line(aes(x = fechas_corregidas, y = tasa_crecimiento_var1, color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas_corregidas, y = tasa_tendencia_var1,  color = "Tendencia"), linewidth = 0.9) +
  scale_color_manual(values = c("Serie original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 1: Tasa de crecimiento anual (%) de la serie original y la tendencia") +
  xlab("Tiempo") + ylab("% de crecimiento anual") +
  theme_minimal()

ggplotly(grafico_crecimiento_var1)

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

tasa_crecimiento_var2 <- (as.numeric(variable2_ts)[13:length(variable2_ts)] /
                            as.numeric(variable2_ts)[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2   <- (as.numeric(tendencia_var2)[13:length(tendencia_var2)] /
                            as.numeric(tendencia_var2)[1:(length(tendencia_var2) - 12)] - 1) * 100

grafico_crecimiento_var2 <- ggplot() +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_line(aes(x = fechas_corregidas, y = tasa_crecimiento_var2, color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas_corregidas, y = tasa_tendencia_var2,  color = "Tendencia"), linewidth = 0.9) +
  scale_color_manual(values = c("Serie original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 2: Tasa de crecimiento anual (%) de la serie original y la tendencia") +
  xlab("Tiempo") + ylab("% de crecimiento anual") +
  theme_minimal()

ggplotly(grafico_crecimiento_var2)

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

tasa_crecimiento_var3 <- (as.numeric(variable3_ts)[13:length(variable3_ts)] /
                            as.numeric(variable3_ts)[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3   <- (as.numeric(tendencia_var3)[13:length(tendencia_var3)] /
                            as.numeric(tendencia_var3)[1:(length(tendencia_var3) - 12)] - 1) * 100

grafico_crecimiento_var3 <- ggplot() +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_line(aes(x = fechas_corregidas, y = tasa_crecimiento_var3, color = "Serie original"), linewidth = 0.5) +
  geom_line(aes(x = fechas_corregidas, y = tasa_tendencia_var3,  color = "Tendencia"), linewidth = 0.9) +
  scale_color_manual(values = c("Serie original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 3: Tasa de crecimiento anual (%) de la serie original y la tendencia") +
  xlab("Tiempo") + ylab("% de crecimiento anual") +
  theme_minimal()

ggplotly(grafico_crecimiento_var3)

Resumen de la tendencia en el último dato disponible (diciembre de 2025)

# Cuántos meses seguidos lleva la tendencia en terreno negativo, contando desde el final
meses_en_negativo <- function(tasa) {
  conteo <- 0
  for (i in length(tasa):1) {
    if (tasa[i] < 0) conteo <- conteo + 1 else break
  }
  conteo
}

resumen_tendencia <- data.frame(
  Variable = c("Cemento Valle", "Licencias Cali", "Cemento Colombia"),
  YoY_dic2025_pct = round(c(tail(tasa_tendencia_var1, 1), tail(tasa_tendencia_var2, 1),
                            tail(tasa_tendencia_var3, 1)), 1),
  Cambio_6m_pp = round(c(tail(tasa_tendencia_var1, 1) - tail(tasa_tendencia_var1, 7)[1],
                         tail(tasa_tendencia_var2, 1) - tail(tasa_tendencia_var2, 7)[1],
                         tail(tasa_tendencia_var3, 1) - tail(tasa_tendencia_var3, 7)[1]), 1),
  Meses_negativos_seguidos = c(meses_en_negativo(tasa_tendencia_var1),
                               meses_en_negativo(tasa_tendencia_var2),
                               meses_en_negativo(tasa_tendencia_var3)),
  Percentil_historico = round(c(100 * mean(tasa_tendencia_var1 <= tail(tasa_tendencia_var1, 1)),
                                100 * mean(tasa_tendencia_var2 <= tail(tasa_tendencia_var2, 1)),
                                100 * mean(tasa_tendencia_var3 <= tail(tasa_tendencia_var3, 1))), 1)
)
resumen_tendencia
Variable YoY_dic2025_pct Cambio_6m_pp Meses_negativos_seguidos Percentil_historico
Cemento Valle -3.7 -2.5 9 11.5
Licencias Cali -77.7 -50.6 8 0.6
Cemento Colombia 9.2 6.0 0 100.0

Interpretación. La tendencia del cemento en el Valle crece -3,7 % anual en diciembre de 2025, después de varios meses consecutivos en terreno negativo, mientras la nacional crece 9,2 %. Las licencias de Cali están en el percentil más bajo de su historia. La divergencia entre el Valle y Colombia sugiere que la causa no sería macroeconómica nacional (tasas de interés, costos, confianza), porque esas afectarían a todo el país por igual, sino específica de la demanda de obra del suroccidente.

Nota metodológica: los filtros de descomposición son asimétricos en el extremo de la muestra, de modo que los últimos 3 a 6 datos de la tendencia se revisarán cuando lleguen nuevas observaciones. La dirección es confiable; el nivel exacto del último punto, menos.

Estacionalidad

meses <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun", "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")

factores <- data.frame(
  Mes = factor(meses, levels = meses),
  Cemento_Valle_pct    = round(as.numeric(estacional_var1)[1:12], 1),
  Licencias_Cali_pct   = round(as.numeric(estacional_var2)[1:12], 1),
  Cemento_Colombia_pct = round(as.numeric(estacional_var3)[1:12], 1)
)
factores
Mes Cemento_Valle_pct Licencias_Cali_pct Cemento_Colombia_pct
Ene -10.7 -14.6 -11.0
Feb -4.7 14.0 -3.3
Mar 1.6 -20.4 5.4
Abr -5.3 9.5 -3.0
May 1.2 -0.1 0.6
Jun -4.0 -11.5 -4.5
Jul 5.7 -10.5 3.4
Ago 6.1 -6.6 3.1
Sep 5.6 -5.6 4.1
Oct 7.0 -17.3 5.6
Nov 4.0 15.5 2.4
Dic -4.6 76.9 -1.5 
grafico_estacional <- ggplot(factores, aes(x = Mes, y = Cemento_Valle_pct)) +
  geom_col(fill = "#1f5f99") +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  ggtitle("Variable 1: Factores estacionales (desviación % frente a la tendencia)") +
  xlab("Mes") + ylab("%") +
  theme_minimal()

ggplotly(grafico_estacional)

Interpretación. El patrón estacional del cemento en el Valle es nítido y estable: el primer trimestre es el más débil del año (enero -10,7 %, febrero -4,7 %) y el segundo semestre el más fuerte (julio a octubre entre +5 % y +7 %). Esto es coherente con el calendario de obra en Colombia: enero concentra vacaciones colectivas y cierres administrativos, mientras el segundo semestre concentra la ejecución presupuestal de obra pública y el cierre de proyectos privados.

Implicación directa para la planta de Yumbo: el trimestre que se va a pronosticar es, por construcción, el punto bajo del año. Una caída de enero frente a diciembre no es una mala noticia por sí misma; la pregunta correcta es si la caída es mayor que el −10,7 % estacional.

En las licencias de Cali el factor de diciembre es muy grande (+76,9 %), pero la dispersión del componente irregular de esa serie es tan alta que ese patrón debe tomarse con cautela: se trata de una serie de aprobaciones “grumosas”, donde un solo proyecto grande mueve el mes.

Componente irregular

irregulares <- data.frame(
  Fecha = format(data_col$FECHA, "%Y-%m"),
  Cemento_Valle_pct = round(as.numeric(irregular_var1), 1)
)

# Vector con nombre de mes, para citar meses concretos en el texto
irr_mes <- setNames(as.numeric(irregular_var1), format(data_col$FECHA, "%Y-%m"))

head(irregulares[order(-abs(irregulares$Cemento_Valle_pct)), ], 6)
Fecha Cemento_Valle_pct
113 2021-05 -87.9
100 2020-04 -75.5
99 2020-03 -42.7
101 2020-05 -34.0
90 2019-06 20.2
15 2013-03 -15.7 
volatilidad <- data.frame(
  Variable = c("Cemento Valle", "Licencias Cali", "Cemento Colombia"),
  DE_irregular_pct = round(c(sd(irregular_var1), sd(irregular_var2), sd(irregular_var3)), 1),
  DE_sin_choques_pct = round(c(
    sd(irregular_var1[!format(data_col$FECHA, "%Y-%m") %in% c("2020-03","2020-04","2020-05","2021-05")]),
    sd(irregular_var2[!format(data_col$FECHA, "%Y-%m") %in% c("2020-03","2020-04","2020-05","2021-05")]),
    sd(irregular_var3[!format(data_col$FECHA, "%Y-%m") %in% c("2020-03","2020-04","2020-05","2021-05")])), 1)
)
volatilidad
Variable DE_irregular_pct DE_sin_choques_pct
Cemento Valle 10.9 4.9
Licencias Cali 65.8 65.9
Cemento Colombia 8.1 3.8 
grafico_irregular <- ggplot() +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_line(aes(x = fechas, y = as.numeric(irregular_var1)), color = "#1f5f99", linewidth = 0.5) +
  ggtitle("Variable 1: Componente irregular (% de desviación frente a tendencia y estacionalidad)") +
  xlab("Tiempo") + ylab("%") +
  theme_minimal()

ggplotly(grafico_irregular)

Interpretación evento por evento:

Mes Desviación irregular Explicación
2020-03 -42,7 % Aislamiento preventivo obligatorio desde el 25 de marzo de 2020.
2020-04 -75,5 % Mes completo de confinamiento; la construcción se reactiva por fases desde el 27 de abril (Decreto 593 de 2020).
2021-05 -87,9 % Paro Nacional: bloqueos en Cali y en los corredores viales del Valle desde el 28 de abril de 2021. Es el mínimo histórico de la serie.

Sin esos cuatro meses, la desviación típica del componente irregular del Valle baja de 10,9 % a 4,9 %. Esto es relevante para el pronóstico: la serie es, fuera de choques, bastante predecible, y por eso los modelos que tratan explícitamente esos meses como atípicos estiman mejor la estructura ordinaria de la serie.

La variable de licencias tiene un componente irregular con desviación típica de 65,8 %, un orden de magnitud mayor. La hipótesis es que las aprobaciones se concentran en pocos proyectos grandes y en cierres de año; verificarlo exige los registros de las Curadurías Urbanas de Cali.

Integración de señales

Las tres series dicen cosas distintas. Integrarlas consiste en responder tres preguntas: ¿la caída del Valle es regional o nacional?, ¿las licencias anticipan al cemento?, y ¿en qué fase del ciclo está hoy cada variable?

# Se reúnen las tasas de crecimiento anual de las tendencias en un solo data frame
senales <- data.frame(
  Fecha            = fechas_corregidas,
  Cemento_Valle    = tasa_tendencia_var1,
  Licencias_Cali   = tasa_tendencia_var2,
  Cemento_Colombia = tasa_tendencia_var3
)
senales$Brecha <- senales$Cemento_Valle - senales$Cemento_Colombia   # diferencia en puntos porcentuales

tail(cbind(Fecha = format(senales$Fecha, "%Y-%m"), round(senales[, -1], 1)), 6)
Fecha Cemento_Valle Licencias_Cali Cemento_Colombia Brecha
151 2025-07 -1.8 -40.6 4.6 -6.4
152 2025-08 -2.2 -51.6 5.8 -8.0
153 2025-09 -2.6 -60.7 7.0 -9.7
154 2025-10 -3.0 -67.4 7.9 -10.9
155 2025-11 -3.3 -72.9 8.8 -12.1
156 2025-12 -3.7 -77.7 9.2 -12.9

1. ¿Regional o nacional? Brecha Valle - Colombia

grafico_brecha <- ggplot(senales, aes(x = Fecha)) +
  geom_col(aes(y = Brecha), fill = "grey70", alpha = 0.6) +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_line(aes(y = Cemento_Valle,    color = "Cemento Valle"), linewidth = 0.9) +
  geom_line(aes(y = Cemento_Colombia, color = "Cemento Colombia"), linewidth = 0.9) +
  scale_color_manual(values = c("Cemento Valle" = "#1f5f99", "Cemento Colombia" = "#7f7f7f")) +
  ggtitle("Crecimiento anual de la tendencia: Valle vs. Colombia (barras = brecha en pp)") +
  xlab("Tiempo") + ylab("% anual / puntos porcentuales") +
  theme_minimal()

ggplotly(grafico_brecha)
cat("Correlación contemporánea (tendencias YoY) Valle - Colombia:",
    round(cor(senales$Cemento_Valle, senales$Cemento_Colombia), 3), "\n")
## Correlación contemporánea (tendencias YoY) Valle - Colombia: 0.569
cat("Brecha promedio 2013-2019 (pp):",
    round(mean(senales$Brecha[senales$Fecha < as.Date("2020-01-01")]), 1), "\n")
## Brecha promedio 2013-2019 (pp): 1
cat("Brecha en diciembre de 2025 (pp):", round(tail(senales$Brecha, 1), 1), "\n")
## Brecha en diciembre de 2025 (pp): -12.9

Interpretación. Históricamente las dos series se mueven juntas (correlación 0,57), pero al cierre de 2025 la brecha se abre a -12,9 puntos porcentuales en contra del Valle. Un choque nacional —costos, tasas de interés, clima— se reflejaría en ambas series; una brecha de esta magnitud apunta a la demanda de obra del suroccidente. La consecuencia para la empresa es directa: el problema no se resuelve esperando a que “el país repunte”, porque el país ya está repuntando.

2. ¿Las licencias anticipan los despachos? Correlación cruzada

La hipótesis del sector es encadenada: se aprueba la licencia, se prevende, se inicia la obra y solo entonces se compra cemento. Si es cierta, la tasa de crecimiento de las licencias de hoy debe parecerse a la del cemento de dentro de varios meses. Eso se mide con la función de correlación cruzada (CCF).

# ccf(x, y) calcula cor(x[t+k], y[t]); si las licencias ADELANTAN al cemento,
# la correlación es alta para rezagos negativos, por eso se invierte el signo del eje.
ccf_simple <- ccf(senales$Licencias_Cali, senales$Cemento_Valle, lag.max = 18, plot = FALSE)

ccf_df <- data.frame(Adelanto = -ccf_simple$lag[, 1, 1], Correlacion = round(ccf_simple$acf[, 1, 1], 3))
ccf_df <- ccf_df[ccf_df$Adelanto >= 0, ]
banda  <- 2 / sqrt(nrow(senales))   # banda de significancia aproximada al 95 %

head(ccf_df[order(-ccf_df$Correlacion), ], 5)
Adelanto Correlacion
8 11 0.379
7 12 0.374
9 10 0.370
6 13 0.355
10 9 0.349

Las dos series son muy persistentes (cada una se parece mucho a sí misma el mes anterior), y eso infla las correlaciones cruzadas. El procedimiento estándar de Box y Jenkins consiste en preblanquear: se ajusta un modelo AR a la serie líder, se filtran ambas series con ese mismo modelo y se calcula la correlación entre lo que queda. Lo que sobreviva al filtro sí es señal.

# Preblanqueo: se ajusta un AR a las licencias y se filtran AMBAS series con ese mismo filtro
ar_licencias <- ar(senales$Licencias_Cali, aic = TRUE, order.max = 12)
filtro       <- c(1, -ar_licencias$ar)

lic_filtrada <- stats::filter(senales$Licencias_Cali - mean(senales$Licencias_Cali), filtro, sides = 1)
cem_filtrada <- stats::filter(senales$Cemento_Valle  - mean(senales$Cemento_Valle),  filtro, sides = 1)
ok <- complete.cases(lic_filtrada, cem_filtrada)

ccf_pb    <- ccf(as.numeric(lic_filtrada)[ok], as.numeric(cem_filtrada)[ok], lag.max = 18, plot = FALSE)
ccf_pb_df <- data.frame(Adelanto = -ccf_pb$lag[, 1, 1], Correlacion = round(ccf_pb$acf[, 1, 1], 3))
ccf_pb_df <- ccf_pb_df[ccf_pb_df$Adelanto >= 0, ]
banda_pb  <- 2 / sqrt(sum(ok))

cat("Orden AR usado en el preblanqueo:", ar_licencias$order,
    "| banda 95 %: +/-", round(banda_pb, 3), "\n")
## Orden AR usado en el preblanqueo: 4 | banda 95 %: +/- 0.162
head(ccf_pb_df[order(-ccf_pb_df$Correlacion), ], 5)
Adelanto Correlacion
7 12 0.239
5 14 0.228
9 10 0.227
11 8 0.227
6 13 0.211 
ccf_grafico <- rbind(
  data.frame(ccf_df,    Tipo = "CCF simple",        Banda = banda),
  data.frame(ccf_pb_df, Tipo = "CCF preblanqueada", Banda = banda_pb)
)

p_ccf <- ggplot(ccf_grafico, aes(x = Adelanto, y = Correlacion)) +
  geom_hline(yintercept = 0, linewidth = 0.3) +
  geom_hline(aes(yintercept = Banda),  linetype = "dashed", color = "#c0392b") +
  geom_hline(aes(yintercept = -Banda), linetype = "dashed", color = "#c0392b") +
  geom_segment(aes(xend = Adelanto, yend = 0), color = "#d9731a", linewidth = 1.1) +
  facet_wrap(~ Tipo, ncol = 2) +
  ggtitle("¿Cuántos meses anticipan las licencias de Cali a los despachos del Valle?") +
  xlab("Adelanto de las licencias, k (meses)") + ylab("Correlación") +
  theme_minimal()

ggplotly(p_ccf)
# El adelanto óptimo se elige con la CCF preblanqueada, que es la que permite probar significancia
k_optimo <- ccf_pb_df$Adelanto[which.max(ccf_pb_df$Correlacion)]

cat("Adelanto óptimo (meses):", k_optimo,
    "| correlación preblanqueada:", round(max(ccf_pb_df$Correlacion), 3),
    "| banda 95 %: +/-", round(banda_pb, 3), "\n")
## Adelanto óptimo (meses): 12 | correlación preblanqueada: 0.239 | banda 95 %: +/- 0.162
cat("Adelantos significativos al 95 %:",
    paste(ccf_pb_df$Adelanto[ccf_pb_df$Correlacion > banda_pb], collapse = ", "), "meses\n")
## Adelantos significativos al 95 %: 15, 14, 13, 12, 11, 10, 9, 8, 7, 6 meses

3. Cuantificación: regresión del cemento sobre las licencias rezagadas

# Se alinean las series: cemento en t contra licencias en t - k_optimo
n_s <- nrow(senales)
datos_lider <- data.frame(
  cemento   = senales$Cemento_Valle[(k_optimo + 1):n_s],
  licencias = senales$Licencias_Cali[1:(n_s - k_optimo)]
)

regresion_lider <- lm(cemento ~ licencias, data = datos_lider)

# Errores estándar robustos a autocorrelación y heterocedasticidad (Newey-West):
# necesarios aquí porque los residuos de una regresión entre dos series de tiempo
# están autocorrelacionados y los errores clásicos serían demasiado optimistas.
coeftest(regresion_lider, vcov = NeweyWest(regresion_lider, lag = 12, prewhite = FALSE))
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.533622   1.137023  1.3488 0.179547   
## licencias   0.112200   0.040414  2.7763 0.006241 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("R2:", round(summary(regresion_lider)$r.squared, 3),
    "| observaciones:", nrow(datos_lider), "| adelanto:", k_optimo, "meses\n")
## R2: 0.219 | observaciones: 144 | adelanto: 12 meses

Interpretación. Un punto porcentual adicional de crecimiento del área licenciada en Cali se asocia con 0,112 puntos porcentuales de crecimiento de los despachos al Valle 12 meses después, y el coeficiente es significativo incluso con errores estándar robustos. El R² de 0,22 dice lo que hay que decir: las licencias informan, pero explican menos de una cuarta parte de la variación. Sirven como alerta temprana, no como pronóstico.

4. Prueba de robustez: causalidad de Granger con datos crudos

Todo lo anterior se calculó sobre tendencias extraídas, que son series suavizadas. Una prueba honesta consiste en repetir la pregunta sobre los datos originales, sin suavizar.

# Variación anual de las series ORIGINALES (diferencia logarítmica de 12 meses)
var_anual_cemento   <- diff(log(variable1_ts), 12) * 100
var_anual_licencias <- diff(log(variable2_ts), 12) * 100

grangertest(var_anual_cemento ~ var_anual_licencias, order = 6)
Res.Df Df F Pr(>F)
137 NA NA NA
143 -6 0.4209815 0.8640557

Resultado nulo, reportado como tal. Con los datos crudos no se rechaza la hipótesis nula: las licencias no “causan en el sentido de Granger” a los despachos. La lectura correcta no es que la relación no exista, sino que la señal queda enterrada bajo el ruido mensual de las licencias (desviación típica del componente irregular de 65,8 %). Solo emerge cuando se suaviza. Esto limita el uso del indicador: sirve para leer la dirección de las tendencias, no para pronosticar meses concretos.

5. Semáforo de señales y proyección indicativa

fase <- function(nivel, cambio3) {
  if (nivel >= 0 && cambio3 >= 0) "Expansión"
  else if (nivel >= 0 && cambio3 <  0) "Desaceleración"
  else if (nivel <  0 && cambio3 <  0) "Contracción"
  else "Recuperación"
}

semaforo <- data.frame(
  Variable = c("Cemento Valle", "Licencias Cali", "Cemento Colombia"),
  YoY_dic2025 = round(c(tail(senales$Cemento_Valle, 1), tail(senales$Licencias_Cali, 1),
                        tail(senales$Cemento_Colombia, 1)), 1),
  Cambio_3m_pp = round(c(tail(senales$Cemento_Valle, 1)    - tail(senales$Cemento_Valle, 4)[1],
                         tail(senales$Licencias_Cali, 1)   - tail(senales$Licencias_Cali, 4)[1],
                         tail(senales$Cemento_Colombia, 1) - tail(senales$Cemento_Colombia, 4)[1]), 1)
)
semaforo$Fase <- NA
for (i in 1:nrow(semaforo)) {
  semaforo$Fase[i] <- fase(semaforo$YoY_dic2025[i], semaforo$Cambio_3m_pp[i])
}
semaforo$Senal <- ifelse(semaforo$Fase == "Expansión", "Verde",
                  ifelse(semaforo$Fase == "Contracción", "Rojo", "Amarillo"))
semaforo
Variable YoY_dic2025 Cambio_3m_pp Fase Senal
Cemento Valle -3.7 -1.1 Contracción Rojo
Licencias Cali -77.7 -17.0 Contracción Rojo
Cemento Colombia 9.2 2.2 Expansión Verde 
# Qué dicen HOY las licencias sobre los próximos meses de despachos
proyeccion_lider <- data.frame(
  Mes = format(seq(as.Date("2026-01-01"), by = "month", length.out = k_optimo), "%Y-%m"),
  Licencias_rezagadas_YoY = round(tail(senales$Licencias_Cali, k_optimo), 1)
)
proyeccion_lider$Cemento_YoY_proyectado <-
  round(coef(regresion_lider)[1] + coef(regresion_lider)[2] * proyeccion_lider$Licencias_rezagadas_YoY, 1)

proyeccion_lider
Mes Licencias_rezagadas_YoY Cemento_YoY_proyectado
2026-01 19.2 3.7
2026-02 19.4 3.7
2026-03 19.6 3.7
2026-04 3.5 1.9
2026-05 -10.4 0.4
2026-06 -27.0 -1.5
2026-07 -40.6 -3.0
2026-08 -51.6 -4.3
2026-09 -60.7 -5.3
2026-10 -67.4 -6.0
2026-11 -72.9 -6.6
2026-12 -77.7 -7.2

Advertencia sobre esta proyección. Es indicativa, no un pronóstico: proviene de una regresión con R² de 0,22 cuya contraparte con datos crudos (Granger) no resultó significativa. Se incluye porque la dirección que señala es información útil para la decisión, aunque el nivel no sea confiable. El pronóstico formal es el de la sección siguiente.

Modelo ARIMA

Un modelo ARIMA (Autoregressive Integrated Moving Average) utiliza el pasado de la propia serie para estimar su futuro.

  1. AR (autorregresivo) → usa valores pasados de la serie.
  2. I (integrado) → diferencia la serie para hacerla estacionaria.
  3. MA (media móvil) → usa los errores de pronóstico pasados.

Cuando la serie repite un patrón cada 12 meses se añade una parte estacional y el modelo se llama SARIMA(p,d,q)(P,D,Q)[12].

Metodología Box-Jenkins

1️⃣ Identificación → 2️⃣ Estimación → 3️⃣ Validación → 4️⃣ Pronóstico.

División en conjunto de entrenamiento y prueba

Se entrena con enero 2012 – septiembre 2025 y se reservan octubre, noviembre y diciembre de 2025 como conjunto de prueba. El horizonte de prueba (3 meses) es el mismo del pronóstico final (enero–marzo 2026), de modo que el error medido corresponde al ejercicio que realmente interesa.

train_ts <- window(variable1_ts, end = c(2025, 9))     # entrenamiento: ene-2012 a sep-2025
test_ts  <- window(variable1_ts, start = c(2025, 10))  # prueba: oct, nov y dic de 2025

cat("Entrenamiento:", length(train_ts), "observaciones |",
    "Prueba:", length(test_ts), "observaciones\n")
## Entrenamiento: 165 observaciones | Prueba: 3 observaciones

Paso 1: Identificación del modelo

Estacionariedad

Una serie es estacionaria si su media y su variabilidad no cambian con el tiempo. Aquí se usan dos pruebas, porque tienen hipótesis nulas opuestas y juntas son mucho más informativas que una sola:

  • ADF (Dickey-Fuller aumentado): H₀ = la serie no es estacionaria. Un p-valor bajo (< 0,05) indica estacionariedad.
  • KPSS: H₀ = la serie sí es estacionaria. Un p-valor bajo (< 0,05) indica que no lo es.
pruebas_raiz <- function(x, nombre) {
  adf  <- suppressWarnings(adf.test(x))
  kpss <- suppressWarnings(kpss.test(x))
  data.frame(Serie = nombre,
             ADF_estadistico  = round(as.numeric(adf$statistic), 3),
             ADF_p            = round(adf$p.value, 3),
             KPSS_estadistico = round(as.numeric(kpss$statistic), 3),
             KPSS_p           = round(kpss$p.value, 3))
}

rbind(
  pruebas_raiz(train_ts,            "Nivel: CEM_V"),
  pruebas_raiz(log(train_ts),       "log(CEM_V)"),
  pruebas_raiz(diff(train_ts),      "Una diferencia"),
  pruebas_raiz(diff(log(train_ts)), "Una diferencia en logaritmo")
)
Serie ADF_estadistico ADF_p KPSS_estadistico KPSS_p
Nivel: CEM_V -5.525 0.01 1.590 0.01
log(CEM_V) -5.737 0.01 1.005 0.01
Una diferencia -7.358 0.01 0.019 0.10
Una diferencia en logaritmo -7.610 0.01 0.017 0.10 
cat("Diferencias regulares sugeridas (ndiffs):", ndiffs(log(train_ts)),
    "| diferencias estacionales sugeridas (nsdiffs):", nsdiffs(log(train_ts)), "\n")
## Diferencias regulares sugeridas (ndiffs): 1 | diferencias estacionales sugeridas (nsdiffs): 0

Lectura de las pruebas. Las dos pruebas se contradicen en el nivel de la serie: la ADF rechaza la raíz unitaria, pero la KPSS también rechaza la estacionariedad. Esto ocurre con frecuencia en series con estacionalidad marcada y atípicos muy grandes: la caída de mayo de 2021 y su recuperación inmediata generan una reversión brusca hacia la media que la ADF interpreta como estacionariedad. Al aplicar una diferencia, ambas pruebas coinciden, y ndiffs() confirma d = 1. Se trabaja entonces con d = 1.

Diferenciación de la variable 1

train_diff     <- diff(train_ts, differences = 1)        # serie diferenciada en niveles
train_diff_log <- diff(log(train_ts), differences = 1)   # serie diferenciada en logaritmo
p2 <- ggplot(data.frame(Tiempo = as.numeric(time(train_ts)), Valor = as.numeric(train_ts)),
             aes(x = Tiempo, y = Valor)) +
  geom_line(color = "blue") +
  ggtitle("Variable 1: Serie original") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal()

ggplotly(p2)
p3 <- ggplot(data.frame(Tiempo = as.numeric(time(train_ts))[-1], Valor = as.numeric(train_diff)),
             aes(x = Tiempo, y = Valor)) +
  geom_line(color = "red") +
  ggtitle("Variable 1: Serie estacionaria (una diferenciación)") +
  xlab("Tiempo") + ylab("Cambio mensual en toneladas") +
  theme_minimal()

ggplotly(p3)

¿Por qué además el logaritmo? Cuando la amplitud de las oscilaciones crece con el nivel de la serie, el logaritmo estabiliza la varianza y convierte los cambios en variaciones porcentuales. En esta serie la diferencia se ve en 2020–2021: sin logaritmo, el choque del Paro Nacional pesa mucho más que cualquier otro movimiento y arrastra la estimación de los parámetros. Por eso los modelos que tratan los atípicos se estiman en logaritmo (lambda = 0).

Identificación manual de p y q con los correlogramas

  • La ACF ayuda a determinar q (parte de media móvil).
  • La PACF ayuda a determinar p (parte autorregresiva).
  • Si una barra sobrepasa las líneas azules, la autocorrelación es significativa al 95 %.
acf_plot  <- ggAcf(train_diff_log, lag.max = 24) +
  ggtitle("ACF de la serie diferenciada en logaritmo - determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 24) +
  ggtitle("PACF de la serie diferenciada en logaritmo - determinar p")

ggplotly(acf_plot)
ggplotly(pacf_plot)

Interpretación de los correlogramas.

  • La PACF presenta barras significativas en los rezagos 1, 2 y 3, y una barra grande en el rezago 12: hay memoria de corto plazo y un componente estacional anual. Valores candidatos: p = 1, 2, 3.
  • La ACF corta rápidamente tras los primeros rezagos: q = 1 o 2.
  • Aparece además una barra significativa en el rezago 13, que no corresponde a ningún ciclo económico: es la distancia exacta entre los dos grandes atípicos de la serie (abril de 2020 y mayo de 2021). Todo apunta a que es un artefacto de los choques, y es la razón por la que más adelante se estima un modelo que los trata explícitamente con variables dummy.

Modelo identificado manualmente: ARIMA(1,1,1), el más parsimonioso compatible con los correlogramas.

Paso 2: Estimación de los modelos candidatos

Se estiman cuatro especificaciones y un referente ingenuo. La comparación con el ingenuo estacional (pronosticar el mismo mes del año anterior) es obligatoria: un modelo que no lo supere no aporta nada.

# Variables dummy para los dos choques identificados en el componente irregular
choques <- cbind(
  covid = as.numeric(format(data_col$FECHA, "%Y-%m") %in% c("2020-03", "2020-04", "2020-05")),
  paro  = as.numeric(format(data_col$FECHA, "%Y-%m") %in% c("2021-05", "2021-06"))
)
xreg_train <- choques[1:length(train_ts), ]
xreg_test  <- choques[(length(train_ts) + 1):(length(train_ts) + length(test_ts)), , drop = FALSE]

Modelo 1: ARIMA(1,1,1) identificado manualmente

modelo_manual <- Arima(train_ts, order = c(1, 1, 1))
summary(modelo_manual)
## Series: train_ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.4920  -0.9443
## s.e.  0.0766   0.0267
## 
## sigma^2 = 114634047:  log likelihood = -1754.02
## AIC=3514.04   AICc=3514.19   BIC=3523.34
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1333.988 10608.95 6622.595 -3.328009 12.03992 0.7376544
##                      ACF1
## Training set -0.000517994
# Significancia estadística de los coeficientes
coeftest(modelo_manual)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1  0.492044   0.076613   6.4224 1.341e-10 ***
## ma1 -0.944269   0.026693 -35.3757 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 2: ARIMA automático sin componente estacional

modelo_auto <- auto.arima(train_ts, seasonal = FALSE)
summary(modelo_auto)
## Series: train_ts 
## ARIMA(0,1,4) 
## 
## Coefficients:
##           ma1      ma2      ma3      ma4
##       -0.4380  -0.2882  -0.0575  -0.1026
## s.e.   0.0785   0.0816   0.0904   0.0739
## 
## sigma^2 = 115644762:  log likelihood = -1753.74
## AIC=3517.48   AICc=3517.86   BIC=3532.98
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 1292.417 10589.64 6641.718 -3.358752 12.02748 0.7397844
##                     ACF1
## Training set -0.01033304
coeftest(modelo_auto)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.437997   0.078525 -5.5778 2.435e-08 ***
## ma2 -0.288163   0.081648 -3.5293 0.0004166 ***
## ma3 -0.057461   0.090409 -0.6356 0.5250635    
## ma4 -0.102646   0.073864 -1.3897 0.1646336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 3: SARIMA automático (con componente estacional)

modelo_sarima <- auto.arima(train_ts)
summary(modelo_sarima)
## Series: train_ts 
## ARIMA(0,1,4)(0,0,1)[12] 
## 
## Coefficients:
##           ma1      ma2      ma3      ma4    sma1
##       -0.4860  -0.2585  -0.0466  -0.1046  0.1289
## s.e.   0.0844   0.0852   0.0915   0.0749  0.0846
## 
## sigma^2 = 114615308:  log likelihood = -1752.57
## AIC=3517.14   AICc=3517.67   BIC=3535.74
## 
## Training set error measures:
##                    ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 1247.275 10509.4 6474.496 -3.504071 11.90838 0.7211585
##                      ACF1
## Training set -0.008326012
coeftest(modelo_sarima)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.485975   0.084438 -5.7554 8.643e-09 ***
## ma2  -0.258493   0.085237 -3.0327  0.002424 ** 
## ma3  -0.046578   0.091498 -0.5091  0.610712    
## ma4  -0.104588   0.074907 -1.3962  0.162643    
## sma1  0.128916   0.084557  1.5246  0.127359    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 4: SARIMA en logaritmo con variables dummy de choque

Los dos atípicos de 2020–2021 no son demanda: son confinamiento y bloqueos. Si no se tratan, el modelo los interpreta como variabilidad normal de la serie y ensancha todos los intervalos de predicción. Aquí se incluyen como regresores externos.

modelo_choques <- auto.arima(train_ts, lambda = 0, biasadj = TRUE, xreg = xreg_train)
summary(modelo_choques)
## Series: train_ts 
## Regression with ARIMA(1,1,1)(0,0,2)[12] errors 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1      ma1    sma1    sma2    covid     paro
##       -0.3545  -0.8600  0.1140  0.1470  -0.8555  -1.1421
## s.e.   0.0832   0.0423  0.0859  0.0734   0.0766   0.0892
## 
## sigma^2 = 0.02381:  log likelihood = 75.54
## AIC=-137.07   AICc=-136.36   BIC=-115.38
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE      ACF1
## Training set 579.2267 9522.778 6875.701 -0.9684334 9.246374 0.7658465 0.3880818
coeftest(modelo_choques)
## 
## z test of coefficients:
## 
##        Estimate Std. Error  z value  Pr(>|z|)    
## ar1   -0.354482   0.083188  -4.2612 2.033e-05 ***
## ma1   -0.860046   0.042325 -20.3200 < 2.2e-16 ***
## sma1   0.114034   0.085862   1.3281    0.1841    
## sma2   0.146978   0.073421   2.0019    0.0453 *  
## covid -0.855506   0.076632 -11.1639 < 2.2e-16 ***
## paro  -1.142053   0.089202 -12.8030 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación de los coeficientes de choque. Como el modelo está en logaritmo, cada coeficiente se convierte en efecto porcentual mediante exp(β) − 1: el confinamiento implicó -57,5 % sobre el nivel esperado de esos meses y el Paro Nacional -68,1 %. Ambos son altamente significativos, lo que confirma que se trata de eventos y no de ruido.

Paso 3: Validación de residuos

El objetivo es verificar que los residuos se comporten como ruido blanco: sin patrones, con media cero y sin autocorrelación. Si queda estructura en los residuos, el modelo está dejando información sin usar.

prueba_residuos <- function(modelo, nombre) {
  # arimaorder devuelve (p, d, q, P, D, Q, m); los grados de libertad consumidos
  # son los parámetros estimados: p + q + P + Q
  k <- sum(arimaorder(modelo)[c(1, 3, 4, 6)], na.rm = TRUE)
  data.frame(
    Modelo         = nombre,
    Especificacion = paste0("(", paste(arimaorder(modelo), collapse = ","), ")"),
    AICc           = round(modelo$aicc, 1),
    LjungBox_12_p  = round(Box.test(residuals(modelo), lag = 12, type = "Ljung-Box", fitdf = k)$p.value, 3),
    LjungBox_24_p  = round(Box.test(residuals(modelo), lag = 24, type = "Ljung-Box", fitdf = k)$p.value, 3),
    JarqueBera_p   = round(jarque.bera.test(residuals(modelo))$p.value, 3)
  )
}

rbind(
  prueba_residuos(modelo_manual,  "M1 ARIMA(1,1,1) manual"),
  prueba_residuos(modelo_auto,    "M2 ARIMA automático"),
  prueba_residuos(modelo_sarima,  "M3 SARIMA automático"),
  prueba_residuos(modelo_choques, "M4 SARIMA log + choques")
)
Modelo Especificacion AICc LjungBox_12_p LjungBox_24_p JarqueBera_p
M1 ARIMA(1,1,1) manual (1,1,1) 3514.2 0.083 0.002 0
M2 ARIMA automático (0,1,4) 3517.9 0.051 0.001 0
M3 SARIMA automático (0,1,4,0,0,1,12) 3517.7 0.144 0.003 0
M4 SARIMA log + choques (1,1,1,0,0,2,12) -136.4 0.217 0.004 0

Lectura honesta del diagnóstico. A 12 rezagos ninguno de los cuatro modelos rechaza la hipótesis de ruido blanco (p entre 0,05 y 0,22), pero a 24 rezagos todos la rechazan, y la prueba de Jarque-Bera rechaza la normalidad en los cuatro casos. Es decir: la estructura de corto plazo está bien capturada y lo que queda sin modelar aparece en los rezagos largos. La causa es la misma que se vio en el correlograma: los dos choques de 2020–2021 generan correlación alrededor del rezago 13 y colas gruesas que ningún ARIMA lineal puede absorber. Las consecuencias prácticas son dos: (i) los pronósticos puntuales siguen siendo utilizables, porque la estructura restante es pequeña; (ii) los intervalos de predicción no están bien calibrados, así que la incertidumbre reportada debe leerse como una aproximación y no como una probabilidad exacta.

checkresiduals(modelo_sarima, lag = 24)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,4)(0,0,1)[12]
## Q* = 40.266, df = 19, p-value = 0.003019
## 
## Model df: 5.   Total lags used: 24

Paso 4: Comparación en el conjunto de prueba (oct - dic 2025)

pron_manual  <- forecast(modelo_manual,  h = length(test_ts))
pron_auto    <- forecast(modelo_auto,    h = length(test_ts))
pron_sarima  <- forecast(modelo_sarima,  h = length(test_ts))
pron_choques <- forecast(modelo_choques, h = length(test_ts), xreg = xreg_test)
pron_ingenuo <- snaive(train_ts, h = length(test_ts))

metricas_prueba <- function(pron, nombre) {
  a <- accuracy(pron, test_ts)
  data.frame(Modelo = nombre,
             RMSE = round(a["Test set", "RMSE"]),
             MAE  = round(a["Test set", "MAE"]),
             MAPE = round(a["Test set", "MAPE"], 2))
}

rbind(
  metricas_prueba(pron_manual,  "M1 ARIMA(1,1,1) manual"),
  metricas_prueba(pron_auto,    "M2 ARIMA automático"),
  metricas_prueba(pron_sarima,  "M3 SARIMA automático"),
  metricas_prueba(pron_choques, "M4 SARIMA log + choques"),
  metricas_prueba(pron_ingenuo, "Referente: ingenuo estacional")
)
Modelo RMSE MAE MAPE
M1 ARIMA(1,1,1) manual 6960 5396 6.50
M2 ARIMA automático 6899 4824 5.89
M3 SARIMA automático 6734 4887 5.94
M4 SARIMA log + choques 6322 5804 6.79
Referente: ingenuo estacional 6824 6215 7.25

Cómo leer las métricas. El MAE es el error promedio en toneladas; el RMSE penaliza más los errores grandes; el MAPE expresa el error como porcentaje del valor real. Como referencia de lectura, un MAPE inferior al 10 % corresponde a un modelo preciso, entre 10 % y 20 % a uno aceptable y por encima de 20 % a uno pobre.

tabla_prueba <- data.frame(
  Mes        = format(seq(as.Date("2025-10-01"), by = "month", length.out = length(test_ts)), "%Y-%m"),
  Observado  = round(as.numeric(test_ts)),
  M1_manual  = round(as.numeric(pron_manual$mean)),
  M2_auto    = round(as.numeric(pron_auto$mean)),
  M3_sarima  = round(as.numeric(pron_sarima$mean)),
  M4_choques = round(as.numeric(pron_choques$mean))
)
tabla_prueba
Mes Observado M1_manual M2_auto M3_sarima M4_choques
2025-10 94708 93401 94799 94865 89900
2025-11 88955 92337 91708 92279 92314
2025-12 80316 91814 91945 91495 89562 
p_prueba <- ggplot(tabla_prueba, aes(x = Mes, group = 1)) +
  geom_line(aes(y = Observado, color = "Observado"), linewidth = 1) +
  geom_point(aes(y = Observado, color = "Observado"), size = 2) +
  geom_line(aes(y = M3_sarima, color = "Pronóstico M3 (SARIMA)"), linewidth = 1) +
  geom_point(aes(y = M3_sarima, color = "Pronóstico M3 (SARIMA)"), size = 2) +
  scale_color_manual(values = c("Observado" = "black", "Pronóstico M3 (SARIMA)" = "#1f5f99")) +
  ggtitle("Variable 1: Pronóstico vs. observado en el conjunto de prueba") +
  xlab("Mes") + ylab("Toneladas") +
  theme_minimal()

ggplotly(p_prueba)

Advertencia. El conjunto de prueba tiene tres observaciones. Elegir un modelo con tres datos es frágil: basta un diciembre atípico para invertir el orden. Por eso la selección no se hace aquí, sino con validación cruzada.

Paso 5: Validación cruzada con origen rodante

La validación cruzada de series de tiempo repite el ejercicio anterior muchas veces: se entrena hasta un mes, se pronostican 1, 2 y 3 meses, se guarda el error, se avanza un mes y se repite. Así cada modelo se evalúa sobre decenas de pronósticos y no sobre tres. La función tsCV() del paquete forecast lo hace automáticamente.

# Cada función recibe la serie hasta un origen y devuelve el pronóstico h pasos adelante.
# Se conservan los órdenes ya identificados y se reestiman los coeficientes en cada origen.
orden_auto   <- arimaorder(modelo_auto)
orden_sarima <- arimaorder(modelo_sarima)
orden_choq   <- arimaorder(modelo_choques)

f_manual  <- function(y, h) forecast(Arima(y, order = c(1, 1, 1)), h = h)
f_auto    <- function(y, h) forecast(Arima(y, order = orden_auto[1:3]), h = h)
f_sarima  <- function(y, h) forecast(Arima(y, order = orden_sarima[1:3], seasonal = orden_sarima[4:6]), h = h)
f_choques <- function(y, h, xreg, newxreg) {
  forecast(Arima(y, order = orden_choq[1:3], seasonal = orden_choq[4:6], lambda = 0, xreg = xreg),
           h = h, xreg = newxreg)
}
f_ingenuo <- function(y, h) snaive(y, h = h)

# initial = 120 -> el primer origen usa 10 años de datos; se evalúa sobre 2022-2025
errores <- list(
  "M1 ARIMA(1,1,1) manual"        = tsCV(variable1_ts, f_manual,  h = 3, initial = 120),
  "M2 ARIMA automático"           = tsCV(variable1_ts, f_auto,    h = 3, initial = 120),
  "M3 SARIMA automático"          = tsCV(variable1_ts, f_sarima,  h = 3, initial = 120),
  "M4 SARIMA log + choques"       = tsCV(variable1_ts, f_choques, h = 3, initial = 120, xreg = choques),
  "Referente: ingenuo estacional" = tsCV(variable1_ts, f_ingenuo, h = 3, initial = 120)
)
# Métricas por horizonte de pronóstico
resumen_cv <- data.frame()

for (nm in names(errores)) {
  e <- errores[[nm]]
  for (h in 1:3) {
    # valor observado que corresponde a cada error del horizonte h
    real <- c(as.numeric(variable1_ts)[-(1:h)], rep(NA, h))
    fila <- data.frame(Modelo = nm, h = h,
                       RMSE = round(sqrt(mean(e[, h]^2, na.rm = TRUE))),
                       MAE  = round(mean(abs(e[, h]), na.rm = TRUE)),
                       MAPE = round(100 * mean(abs(e[, h] / real), na.rm = TRUE), 2),
                       n    = sum(!is.na(e[, h])))
    resumen_cv <- rbind(resumen_cv, fila)
  }
}
resumen_cv[order(resumen_cv$h, resumen_cv$RMSE), ]
Modelo h RMSE MAE MAPE n
1 M1 ARIMA(1,1,1) manual 1 7608 6344 6.85 47
4 M2 ARIMA automático 1 7767 6428 6.94 47
7 M3 SARIMA automático 1 8565 6563 7.04 47
10 M4 SARIMA log + choques 1 11783 8556 9.23 47
13 Referente: ingenuo estacional 1 17539 10724 11.14 47
2 M1 ARIMA(1,1,1) manual 2 8977 7236 7.75 46
5 M2 ARIMA automático 2 9235 7369 7.89 46
8 M3 SARIMA automático 2 9814 7497 7.99 46
11 M4 SARIMA log + choques 2 13161 9034 9.66 46
14 Referente: ingenuo estacional 2 17661 10728 11.15 46
3 M1 ARIMA(1,1,1) manual 3 9608 7583 8.14 45
6 M2 ARIMA automático 3 9773 7696 8.27 45
9 M3 SARIMA automático 3 10280 7829 8.37 45
12 M4 SARIMA log + choques 3 13145 9228 9.94 45
15 Referente: ingenuo estacional 3 17739 10663 11.11 45 
# Error global y error restringido a los meses que se van a pronosticar (enero, febrero y marzo)
resumen_q1 <- data.frame()

for (nm in names(errores)) {
  e <- errores[[nm]]

  # marca los errores cuyo mes objetivo cae en enero, febrero o marzo
  sel <- matrix(FALSE, nrow(e), 3)
  for (h in 1:3) {
    mes_objetivo <- (1:nrow(e) + h - 1) %% 12 + 1
    sel[, h] <- mes_objetivo %in% c(1, 2, 3)
  }

  fila <- data.frame(Modelo = nm,
                     RMSE_global  = round(sqrt(mean(e^2, na.rm = TRUE))),
                     RMSE_ene_mar = round(sqrt(mean(e[sel]^2, na.rm = TRUE))))
  resumen_q1 <- rbind(resumen_q1, fila)
}
resumen_q1[order(resumen_q1$RMSE_global), ]
Modelo RMSE_global RMSE_ene_mar
M1 ARIMA(1,1,1) manual 8757 10467
M2 ARIMA automático 8951 10526
M3 SARIMA automático 9568 9662
M4 SARIMA log + choques 12703 9915
Referente: ingenuo estacional 17645 9339 
p_cv <- ggplot(resumen_cv, aes(x = factor(h), y = MAPE, fill = Modelo)) +
  geom_col(position = position_dodge(width = 0.85), width = 0.8) +
  ggtitle("Validación cruzada con origen rodante (objetivos 2022-2025)") +
  xlab("Horizonte de pronóstico (meses)") + ylab("MAPE (%)") +
  theme_minimal() +
  theme(legend.position = "bottom")

ggplotly(p_cv)

Paso 6: Selección del modelo

El criterio se declara antes de mirar los resultados: se elige el modelo con menor error de validación cruzada en los meses que efectivamente se van a pronosticar (enero a marzo), porque esa es la pérdida relevante para la decisión. El error global y el conjunto de prueba se usan como controles.

Criterio Resultado
RMSE de validación cruzada, meses enero–marzo M3 SARIMA automático (9.662 t); luego M4 (9.915) y M1 (10.467)
RMSE de validación cruzada, global (h = 1 a 3) M1 (8.757 t), M2 (8.951), M3 (9.568), M4 (12.703)
RMSE en el conjunto de prueba (oct–dic 2025) M4 (6.322 t), M3 (6.734), M2 (6.899), M1 (6.960)
AICc entre los modelos estimados en niveles M1 (3.514,2), M3 (3.517,7), M2 (3.517,9)

Modelo seleccionado: M3, el SARIMA automático. Gana en el criterio declarado y queda segundo en el conjunto de prueba. Conviene ser transparente sobre lo que esta tabla muestra y lo que no:

  • Ningún modelo domina a los demás. M1, el identificado manualmente, es el mejor en el error global y en AICc; M3 es el mejor en los meses objetivo; M4 es el mejor en el conjunto de prueba. Las diferencias entre ellos —entre 8.757 y 9.568 toneladas de RMSE global— son del orden del 9 %, dentro del ruido de una validación con unos 45 orígenes. La conclusión sustantiva no depende de cuál se elija, como se verifica más adelante.
  • El ingenuo estacional tiene el menor RMSE de enero a marzo (9.339 t), aunque es, con diferencia, el peor en el conjunto del año (17.645 t frente a 8.757–9.568 de los ARIMA). La explicación más plausible es que en el primer trimestre la serie se mueva poco de un año a otro, y repetir el mismo mes del año anterior es difícil de superar. Es un aviso metodológico de que el valor agregado del modelo en ese trimestre concreto es menor que en el resto del año.
  • M4, con dummies de choque, es el peor en validación cruzada pese a ser el mejor en el conjunto de prueba y el que mejor interpreta económicamente los atípicos. En esta serie, el tratamiento explícito de los choques mejora la estimación estructural pero, combinado con la transformación logarítmica, deteriora la precisión predictiva fuera de muestra. Es un buen ejemplo de que el mejor modelo para explicar no siempre es el mejor modelo para pronosticar.

Paso 7: Pronóstico enero - marzo 2026

El modelo seleccionado se reestima con toda la muestra (enero 2012 – diciembre 2025): para pronosticar el futuro no tiene sentido descartar los últimos tres meses de información.

modelo_final <- Arima(variable1_ts, order = orden_sarima[1:3], seasonal = orden_sarima[4:6])
summary(modelo_final)
## Series: variable1_ts 
## ARIMA(0,1,4)(0,0,1)[12] 
## 
## Coefficients:
##           ma1      ma2      ma3      ma4    sma1
##       -0.4846  -0.2593  -0.0493  -0.1018  0.1303
## s.e.   0.0837   0.0845   0.0888   0.0734  0.0836
## 
## sigma^2 = 113111511:  log likelihood = -1783.56
## AIC=3579.12   AICc=3579.64   BIC=3597.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 1139.879 10443.75 6434.306 -3.543405 11.7875 0.7209484
##                      ACF1
## Training set -0.007110358
coeftest(modelo_final)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1  -0.484599   0.083712 -5.7889 7.085e-09 ***
## ma2  -0.259280   0.084503 -3.0683  0.002153 ** 
## ma3  -0.049261   0.088773 -0.5549  0.578952    
## ma4  -0.101776   0.073434 -1.3860  0.165761    
## sma1  0.130286   0.083590  1.5586  0.119082    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
H <- 3   # enero, febrero y marzo de 2026
pronostico_final <- forecast(modelo_final, h = H, level = c(80, 95))

base_2025      <- as.numeric(window(variable1_ts, start = c(2025, 1), end = c(2025, 3)))  # ene-mar 2025
total_2025     <- sum(base_2025)
fechas_futuras <- seq(as.Date("2026-01-01"), by = "month", length.out = H)

tabla_pronostico <- data.frame(
  Mes            = format(fechas_futuras, "%Y-%m"),
  Pronostico_t   = round(as.numeric(pronostico_final$mean)),
  LI_80          = round(as.numeric(pronostico_final$lower[, 1])),
  LS_80          = round(as.numeric(pronostico_final$upper[, 1])),
  LI_95          = round(as.numeric(pronostico_final$lower[, 2])),
  LS_95          = round(as.numeric(pronostico_final$upper[, 2])),
  Mismo_mes_2025 = round(base_2025),
  Var_anual_pct  = round(100 * (as.numeric(pronostico_final$mean) / base_2025 - 1), 1)
)
tabla_pronostico
Mes Pronostico_t LI_80 LS_80 LI_95 LS_95 Mismo_mes_2025 Var_anual_pct
2026-01 83912 70282 97542 63067 104757 77510 8.3
2026-02 87097 71763 102431 63646 110548 83075 4.8
2026-03 88942 73216 104668 64892 112993 90124 -1.3 
historico <- data.frame(Fecha = fechas, Valor = as.numeric(variable1_ts))
historico <- historico[historico$Fecha >= as.Date("2022-01-01"), ]

futuro <- data.frame(
  Fecha = fechas_futuras,
  Media = as.numeric(pronostico_final$mean),
  LI80  = as.numeric(pronostico_final$lower[, 1]), LS80 = as.numeric(pronostico_final$upper[, 1]),
  LI95  = as.numeric(pronostico_final$lower[, 2]), LS95 = as.numeric(pronostico_final$upper[, 2])
)

p_pronostico <- ggplot() +
  geom_ribbon(data = futuro, aes(x = Fecha, ymin = LI95, ymax = LS95), fill = "#1f5f99", alpha = 0.15) +
  geom_ribbon(data = futuro, aes(x = Fecha, ymin = LI80, ymax = LS80), fill = "#1f5f99", alpha = 0.30) +
  geom_line(data = historico, aes(x = Fecha, y = Valor), color = "grey30", linewidth = 0.6) +
  geom_line(data = futuro, aes(x = Fecha, y = Media), color = "#1f5f99", linewidth = 1) +
  geom_point(data = futuro, aes(x = Fecha, y = Media), color = "#1f5f99", size = 2) +
  ggtitle("Pronóstico de despachos de cemento al Valle, enero - marzo 2026") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal()

ggplotly(p_pronostico)

Robustez: ¿cambiaría la conclusión con otro de los modelos candidatos?

comparacion_final <- data.frame(
  Modelo = c("M1 ARIMA(1,1,1) manual", "M2 ARIMA automático", "M3 SARIMA automático (seleccionado)"),
  Trimestre_t = c(
    sum(forecast(Arima(variable1_ts, order = c(1, 1, 1)), h = H)$mean),
    sum(forecast(Arima(variable1_ts, order = orden_auto[1:3]), h = H)$mean),
    sum(as.numeric(pronostico_final$mean)))
)
comparacion_final$Var_anual_pct <- round(100 * (comparacion_final$Trimestre_t / total_2025 - 1), 1)
comparacion_final$Trimestre_t   <- round(comparacion_final$Trimestre_t)
comparacion_final
Modelo Trimestre_t Var_anual_pct
M1 ARIMA(1,1,1) manual 261586 4.3
M2 ARIMA automático 261560 4.3
M3 SARIMA automático (seleccionado) 259951 3.7

Los tres modelos proyectan un primer trimestre por encima del de 2025, con diferencias de pocos puntos porcentuales entre ellos. La elección del modelo no cambia el mensaje: todos extrapolan continuidad, y esa es precisamente la limitación que la validación ex post pondrá en evidencia.

Paso 8: Tratamiento de la incertidumbre

Un pronóstico puntual es la peor manera de comunicar un resultado a una gerencia, porque sugiere una precisión que no existe. Dos herramientas complementarias:

  1. Los intervalos de predicción del modelo (tabla anterior).
  2. Una simulación de trayectorias futuras, que permite responder preguntas de negocio del tipo “¿cuál es la probabilidad de que el trimestre sea peor que el del año pasado?”.
set.seed(2026)
simulaciones <- replicate(3000, sum(simulate(modelo_final, nsim = H, future = TRUE)))

resumen_trimestre <- data.frame(
  Concepto = c("I trimestre 2025 (observado)", "I trimestre 2026 (pronóstico medio)",
               "Límite inferior 80 %", "Límite superior 80 %",
               "Límite inferior 95 %", "Límite superior 95 %"),
  Toneladas = round(c(total_2025, mean(simulaciones),
                      quantile(simulaciones, 0.10), quantile(simulaciones, 0.90),
                      quantile(simulaciones, 0.025), quantile(simulaciones, 0.975)))
)
resumen_trimestre
Concepto Toneladas
I trimestre 2025 (observado) 250710
I trimestre 2026 (pronóstico medio) 260133
Límite inferior 80 % 226072
Límite superior 80 % 293862
Límite inferior 95 % 206318
Límite superior 95 % 311844 
prob_caida <- mean(simulaciones < total_2025)

cat("Variación anual esperada del I trimestre de 2026:",
    round(100 * (mean(simulaciones) / total_2025 - 1), 1), "%\n")
## Variación anual esperada del I trimestre de 2026: 3.8 %
cat("Probabilidad de que el trimestre sea inferior al de 2025:", round(100 * prob_caida), "%\n")
## Probabilidad de que el trimestre sea inferior al de 2025: 36 %

Interpretación económica del pronóstico. El modelo proyecta un primer trimestre de 259.951 toneladas, un 3,7 % frente al mismo trimestre de 2025, con una probabilidad de caída anual de 36 %. Conviene ser explícito sobre de dónde viene ese resultado: un ARIMA extrapola la inercia reciente de la propia serie y no sabe nada de licencias, tasas de interés ni iniciaciones de vivienda. Como el último trimestre de 2025 fue relativamente alto, el modelo proyecta continuidad. La señal del indicador líder apunta exactamente en la dirección contraria. Ese conflicto entre las dos señales es, en sí mismo, un hallazgo de gestión, y la sección siguiente lo resuelve con datos publicados después del corte de la base.

Validación ex post con datos publicados después del corte

La base de este análisis termina en diciembre de 2025. Al cierre de este informe, en septiembre de 2026, el DANE ya había publicado los despachos reales del primer trimestre de 2026, lo que permite contrastar el pronóstico con el dato observado.

Los boletines de Estadísticas de Cemento Gris reportan, para los despachos con destino al Valle del Cauca (DANE, 2026a, 2026b): enero de 2026 −0,9 % anual, marzo de 2026 −5,6 % anual y el acumulado enero–marzo −2,6 %. Febrero se deriva por diferencia con el acumulado trimestral.

# Variaciones anuales publicadas por el DANE para el Valle del Cauca
var_dane <- c(enero = -0.9, marzo = -5.6, trimestre = -2.6)

real_ene  <- base_2025[1] * (1 + var_dane["enero"] / 100)
real_mar  <- base_2025[3] * (1 + var_dane["marzo"] / 100)
real_trim <- total_2025   * (1 + var_dane["trimestre"] / 100)
real_feb  <- real_trim - real_ene - real_mar     # se deriva por diferencia

reales <- as.numeric(c(real_ene, real_feb, real_mar))

ex_post <- data.frame(
  Mes          = tabla_pronostico$Mes,
  Pronostico_t = tabla_pronostico$Pronostico_t,
  Real_DANE_t  = round(reales),
  Error_t      = round(reales - tabla_pronostico$Pronostico_t),
  Error_pct    = round(100 * (reales - tabla_pronostico$Pronostico_t) / reales, 1),
  Dentro_IC80  = reales >= tabla_pronostico$LI_80 & reales <= tabla_pronostico$LS_80,
  Dentro_IC95  = reales >= tabla_pronostico$LI_95 & reales <= tabla_pronostico$LS_95
)
ex_post
Mes Pronostico_t Real_DANE_t Error_t Error_pct Dentro_IC80 Dentro_IC95
2026-01 83912 76813 -7099 -9.2 TRUE TRUE
2026-02 87097 82301 -4796 -5.8 TRUE TRUE
2026-03 88942 85078 -3864 -4.5 TRUE TRUE 
cat("Trimestre pronosticado:", round(sum(tabla_pronostico$Pronostico_t)), "t\n")
## Trimestre pronosticado: 259951 t
cat("Trimestre real (DANE): ", round(real_trim), "t\n")
## Trimestre real (DANE):  244191 t
cat("Error del trimestre:   ", round(100 * (real_trim - sum(tabla_pronostico$Pronostico_t)) / real_trim, 1), "%\n")
## Error del trimestre:    -6.5 %
cat("¿El trimestre real cae dentro del intervalo simulado al 80 %?:",
    real_trim >= quantile(simulaciones, 0.10) & real_trim <= quantile(simulaciones, 0.90), "\n")
## ¿El trimestre real cae dentro del intervalo simulado al 80 %?: TRUE
reales_df <- data.frame(Fecha = fechas_futuras, Valor = reales)

p_expost <- ggplot() +
  geom_ribbon(data = futuro, aes(x = Fecha, ymin = LI95, ymax = LS95), fill = "#1f5f99", alpha = 0.15) +
  geom_ribbon(data = futuro, aes(x = Fecha, ymin = LI80, ymax = LS80), fill = "#1f5f99", alpha = 0.30) +
  geom_line(data = historico, aes(x = Fecha, y = Valor), color = "grey30", linewidth = 0.6) +
  geom_line(data = futuro, aes(x = Fecha, y = Media), color = "#1f5f99", linewidth = 1) +
  geom_point(data = futuro, aes(x = Fecha, y = Media), color = "#1f5f99", size = 2) +
  geom_point(data = reales_df, aes(x = Fecha, y = Valor), color = "#c0392b", size = 3, shape = 4, stroke = 1.3) +
  ggtitle("Pronóstico (azul) vs. dato real publicado por el DANE (X roja)") +
  xlab("Tiempo") + ylab("Toneladas") +
  theme_minimal()

ggplotly(p_expost)

Qué salió bien y qué salió mal.

  • El modelo pronosticó un trimestre de 259.951 t (+3,7 % anual) y el dato real publicado por el DANE (2026a, 2026b) fue de 244.191 t (−2,6 %). El modelo sobreestimó la demanda en 15.760 toneladas, un 6,5 %.
  • Los datos reales caen dentro del intervalo de predicción al 80 % en los tres meses, de modo que el error está dentro de lo que el propio modelo anunciaba como posible. El intervalo cumplió su función; la media no.
  • El error tiene una dirección clara y no es aleatorio: un ARIMA solo extrapola la inercia de la serie, y esa inercia venía de un cuarto trimestre de 2025 relativamente bueno. La información que anticipaba la caída —licencias en mínimos históricos— estaba fuera del modelo. Es exactamente el tipo de error que un indicador líder está destinado a corregir.
  • Advertencia de comparabilidad: el dato real se reconstruyó aplicando las variaciones publicadas a la base de 2025 de este archivo, y el DANE advirtió en su boletín de enero de 2026 que revisó cifras de 2024–2025 por información extemporánea. La comparación de tasas es sólida; la de niveles mezcla dos versiones de la serie. Para cerrarla habría que descargar el anexo del boletín de marzo de 2026 y comparar niveles del mismo vintage.

¿Y el indicador líder acertó? Conviene ser igual de exigente con él que con el modelo. La proyección indicativa apuntaba a 3,7 % para enero de 2026: también falló en el signo, y por la misma razón que el ARIMA, porque con un adelanto de 12 meses el dato que alimenta ese pronóstico es el pico de licencias de finales de 2024, que nunca se convirtió en iniciaciones de obra. Licencias aprobadas no son obras iniciadas.

Donde sí acertó fue en el mediano plazo. La misma proyección pasa a terreno negativo desde mediados de 2026 y llega a -7,2 % al cierre del año. La evidencia posterior lo confirma: según el boletín del DANE de julio de 2026 (2026c), los despachos al Valle cayeron −11,5 % anual y el acumulado enero–julio −4,3 %. La lección metodológica es doble: las licencias avisan la dirección del ciclo con varios meses de anticipación, pero no permiten calibrar el nivel de un mes concreto, tal como ya anticipaban el R² de 0,22 y el resultado nulo de la prueba de Granger.

Una precisión importante sobre la composición de la caída. Según los boletines del DANE (2026a, 2026b, 2026c), en el Valle el canal granel cayó −26,2 % en enero y −19,2 % en marzo de 2026, mientras el empacado creció +9,7 % y +0,1 %. La contracción está concentrada en la obra formal —concreteras, constructoras y contratistas—, no en el consumo de hogares y ferreterías. Esa distinción es la que convierte el diagnóstico en decisión comercial.

De los resultados a la decisión

Riesgos

Riesgo Evidencia Implicación para la operación de Yumbo
Contracción de la obra formal en el Valle Tendencia del Valle en -3,7 % anual frente a 9,2 % nacional; granel −26,2 % en enero de 2026 (DANE, 2026a) La caída es del segmento de mayor volumen por cliente; afecta utilización de planta y logística de granel
Licencias de Cali en mínimos históricos Área licenciada −36,6 % en 2025; tendencia en -77,7 % anual Con un adelanto estimado de 12 meses, la demanda del cierre de 2026 y de comienzos de 2027 ya está comprometida a la baja
Sobrecapacidad del sector Capacidad instalada nacional ~23,6 Mt frente a despachos ~12,5 Mt promedio de cinco años (Portafolio, 2026) Presión a la baja sobre precios e incentivo a exportar excedentes o reducir turnos
Costo y disponibilidad de energía El sector consume cerca de 6,6 veces el promedio manufacturero; riesgo de fenómeno de El Niño en el segundo semestre de 2026 (Portafolio, 2026) Un choque de costos sobre márgenes ya presionados por menor volumen
Consolidación de la competencia Acuerdo de venta de activos de Cemex a Holcim en Colombia por USD 485 millones (Portafolio, 2026) Cambio en la estructura competitiva regional; falta verificar el estado de aprobación regulatoria

Oportunidades

Oportunidad Evidencia Acción que habilita
El canal empacado resiste Empacado +9,7 % en enero de 2026 en el Valle mientras granel cae −26,2 % (DANE, 2026a) Reasignar esfuerzo comercial y logístico hacia ferretería, distribución y remodelación
Estacionalidad predecible y estable Primer trimestre -4,6 % frente a la tendencia; julio a noviembre entre +4 % y +7 % Programar mantenimientos y vacaciones en el valle estacional, y reservar capacidad para julio–octubre
Mercado nacional en crecimiento Tendencia nacional en 9,2 % anual; despachos nacionales +5,7 % en el primer trimestre de 2026 (DANE, 2026b) Desviar excedentes de Yumbo hacia mercados con crecimiento, si el costo logístico lo permite
Indicador líder disponible y público Adelanto de 12 meses, significativo con errores robustos Incorporar las licencias de Cali al proceso de planeación de la demanda

Decisiones concretas recomendadas

  1. Planear producción e inventarios en el borde inferior del intervalo, no en la media. La validación ex post mostró que el error del modelo fue sistemático y hacia arriba. Usar el límite inferior al 80 % del trimestre simulado (226.072 t) como escenario base de producción y la media (260.133 t) como escenario optimista. El dato real del trimestre, 244.191 t, quedó entre ese borde inferior y la media: planear en el borde habría evitado el exceso de inventario sin generar desabastecimiento.
  2. Combinar el modelo con el indicador líder mediante una regla explícita. Si la variación anual del área licenciada en Cali es menor a −25 % durante tres meses consecutivos, ajustar a la baja el pronóstico del ARIMA en el horizonte de 12 meses hacia adelante. Con los datos disponibles a diciembre de 2025 —licencias en -77,7 %—, la regla no habría corregido el primer trimestre de 2026 (el adelanto apunta al cierre de 2026), pero sí habría anticipado el deterioro observado desde mitad de año.
  3. Recomponer la mezcla comercial de granel a empacado. El diferencial de comportamiento entre canales (−26,2 % frente a +9,7 % en enero) indica dónde está la demanda que sí existe: ferreterías, distribuidores y remodelación.
  4. Programar mantenimientos mayores en enero y febrero. Son los meses estacionalmente más débiles (-10,7 % y -4,7 % frente a la tendencia): el costo de oportunidad de parar es el más bajo del año.
  5. Revisar la política de crédito a distribuidores. Si el canal empacado se vuelve el soporte del volumen, la exposición de cartera se atomiza; conviene monitorear días de cartera por canal, lo que requiere datos internos de la empresa.

Indicadores a monitorear, con umbral

Indicador Fuente Frecuencia Umbral de alerta
Área licenciada en Cali, variación anual DANE – ELIC Mensual Por debajo de −25 % durante 3 meses seguidos
Despachos al Valle, granel vs. empacado DANE – ECG Mensual Granel cayendo más de 15 % anual
Participación del Valle en despachos nacionales DANE – ECG Mensual Por debajo de 8,0 %
Brecha de tendencia Valle − Colombia Cálculo propio de este informe Mensual Más negativa que −10 pp
Iniciaciones de vivienda en el Valle DANE – Censo de Edificaciones Trimestral Caída superior a 20 % anual
Tasa de política monetaria y desembolsos hipotecarios Banco de la República Mensual Umbral sin definir: requiere la trayectoria de tasas de 2026

Limitaciones y vacíos de información

  1. Extremo de la muestra. Los últimos tres a seis datos de la tendencia extraída se revisarán cuando lleguen nuevas observaciones, porque los filtros son asimétricos al final de la serie. La dirección es confiable; el nivel exacto del último punto, menos.
  2. Residuos que no son ruido blanco. Los cuatro modelos rechazan la prueba de Ljung-Box a 24 rezagos por efecto de los dos choques de 2020–2021. Los intervalos de predicción deben interpretarse como aproximados.
  3. Conjunto de prueba muy pequeño. Tres observaciones no permiten seleccionar un modelo; por eso la selección se apoyó en la validación cruzada, con aproximadamente 45 orígenes.
  4. Mezcla de versiones de la serie en la validación ex post. El dato real se reconstruyó con tasas publicadas y el DANE revisó cifras de 2024–2025. La comparación de tasas es válida; la de niveles es aproximada.
  5. Empresa. No se verificaron capacidad instalada, volúmenes ni resultados de la planta de Yumbo. Toda la lectura empresarial se apoya en la demanda regional observada, no en información de la compañía.
  6. Pendientes de verificación. Quedan sin confirmar: licencias e iniciaciones en Cali durante 2026; tasa de política monetaria del Banco de la República en 2026 y su transmisión al crédito hipotecario; ajustes a los subsidios de vivienda de interés social; obras de infraestructura en el Valle con demanda de cemento a granel.
  7. Causalidad. Nada de lo presentado identifica relaciones causales. La relación licencias → despachos es una asociación estadística coherente con el encadenamiento productivo del sector, respaldada por una CCF preblanqueada y una regresión con errores robustos, pero contradicha por la prueba de Granger con datos crudos.

Conclusión

La descomposición de las tres series muestra un mercado regional en contracción (-3,7 % anual en la tendencia al cierre de 2025) dentro de un país en expansión (9,2 %), con un patrón estacional estable que sitúa el primer trimestre como el punto bajo del año y con un componente irregular dominado por dos eventos plenamente identificables: el confinamiento de 2020 y el Paro Nacional de 2021.

El modelo seleccionado por validación cruzada, un SARIMA automático, superó con claridad al referente ingenuo, pero sobreestimó el primer trimestre de 2026 en cerca de un 6,5 %, aunque el dato real se mantuvo dentro de su intervalo al 80 %. El indicador líder construido con las licencias de Cali —significativo, con un adelanto estimado de 12 meses, pero con un R² modesto y sin respaldo de la prueba de Granger sobre datos crudos— tampoco acertó el primer trimestre, aunque sí anticipó la dirección del deterioro que se confirmó a mitad de 2026.

La conclusión de gestión es más útil que la estadística: para la planta de Yumbo, un modelo univariado bien construido sirve para dimensionar el rango de la demanda del próximo trimestre, pero no para anticipar un cambio de ciclo. Eso exige mirar hacia arriba en la cadena —licencias e iniciaciones— y hacia adentro de la mezcla de canales, donde la caída del granel y la resistencia del empacado indican con precisión dónde ajustar producción y dónde concentrar el esfuerzo comercial.

Referencias

Box, G. E. P., & Jenkins, G. M. (1976). Time series analysis: Forecasting and control (ed. rev.). Holden-Day.

Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3–73.

Departamento Administrativo Nacional de Estadística. (2026a). Estadísticas de cemento gris (ECG): enero de 2026 [Boletín técnico]. https://www.dane.gov.co/files/operaciones/ECG/bol-ECG-ene2026.pdf

Departamento Administrativo Nacional de Estadística. (2026b). Estadísticas de cemento gris (ECG): marzo de 2026 [Boletín técnico]. https://www.dane.gov.co/files/operaciones/ECG/bol-ECG-mar2026.pdf

Departamento Administrativo Nacional de Estadística. (2026c). Estadísticas de cemento gris (ECG): julio de 2026 [Boletín técnico]. https://www.dane.gov.co/files/operaciones/ECG/bol-ECG-jul2026.pdf

Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and practice (3.ª ed.). OTexts. https://otexts.com/fpp3/

Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(3), 1–22. https://doi.org/10.18637/jss.v027.i03

Newey, W. K., & West, K. D. (1987). A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica, 55(3), 703–708. https://doi.org/10.2307/1913610

Portafolio. (2026, 4 de junio). Bajón en construcción tiene a cementeras con exceso de capacidad e inventarios al límite en Colombia. https://www.portafolio.co/economia/infraestructura/bajon-en-construccion-tiene-a-cementeras-con-exceso-de-capacidad-e-inventarios-al-limite-en-colombia-495527

Presidencia de la República de Colombia. (2020). Decreto 593 de 2020, por el cual se imparten instrucciones en virtud de la emergencia sanitaria generada por la pandemia del coronavirus COVID-19.

sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] plotly_4.12.1   ggplot2_4.0.3   sandwich_3.1-1  lmtest_0.9-40  
## [5] zoo_1.8-14      forecast_9.0.2  tseries_0.10-61 readxl_1.4.5   
## 
## loaded via a namespace (and not attached):
##  [1] tidyr_1.3.1        sass_0.4.9         utf8_1.2.4         generics_0.1.3    
##  [5] lattice_0.22-6     digest_0.6.36      magrittr_2.0.3     evaluate_1.0.5    
##  [9] grid_4.4.1         RColorBrewer_1.1-3 fastmap_1.2.0      cellranger_1.1.0  
## [13] jsonlite_2.0.0     httr_1.4.7         purrr_1.2.2        fansi_1.0.6       
## [17] crosstalk_1.2.1    viridisLite_0.4.2  scales_1.4.0       jquerylib_0.1.4   
## [21] cli_3.6.3          rlang_1.2.0        withr_3.0.2        cachem_1.1.0      
## [25] yaml_2.3.10        otel_0.2.0         tools_4.4.1        parallel_4.4.1    
## [29] dplyr_1.1.4        colorspace_2.1-1   curl_7.1.0         vctrs_0.6.5       
## [33] R6_2.6.1           lifecycle_1.0.4    htmlwidgets_1.6.4  urca_1.3-4        
## [37] pkgconfig_2.0.3    bslib_0.8.0        pillar_1.9.0       gtable_0.3.6      
## [41] data.table_1.16.0  glue_1.7.0         quantmod_0.4.28    Rcpp_1.0.13       
## [45] xfun_0.56          tibble_3.2.1       tidyselect_1.2.1   rstudioapi_0.17.1 
## [49] knitr_1.51         farver_2.1.2       nlme_3.1-164       htmltools_0.5.8.1 
## [53] labeling_0.4.3     rmarkdown_2.30     xts_0.14.1         timeDate_4052.112 
## [57] fracdiff_1.5-3     compiler_4.4.1     quadprog_1.5-8     S7_0.2.1          
## [61] TTR_0.24.4