Introducción

Para realizar la extracción de señales y la aplicación de los modelos ARIMA y SARIMA del caso 2, se selecciona el sector de construcción, tomando como referencia las siguientes variables:

  1. Despachos de cemento (DECEM)
  2. Producción de cemento (PNCEM)
  3. Licencias de construcción en Colombia (LICC)

Variable a pronosticar: Despachos de cemento (DECEM)

Justificación de la selección

Se selecciona el sector de la construcción y como empresa de referencia a Cementos San Marcos, debido a su relación directa con la producción y comercialización de cemento. Las variables seleccionadas permiten analizar la demanda, la oferta y la actividad constructiva del sector. Se espera identificar relaciones entre su comportamiento para anticipar cambios en la demanda y apoyar decisiones de producción, inventarios y comercialización.

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

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

Variable 1

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

Variable 2

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

Variable 3

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

Extracción 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("V1 DECEM: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal() + 
  scale_y_continuous(labels = scales::comma)

ggplotly(grafico_serie)

Interpretación serie original

La serie original de despachos de cemento presenta un comportamiento fluctuante a lo largo del período analizado, con variaciones importantes entre los diferentes meses. Sin embargo, se observa una tendencia general creciente en el nivel de despachos. Se destaca una caída excepcional alrededor del año 2020, donde los despachos disminuyen significativamente respecto a su comportamiento habitual. Esta caída coincide temporalmente con el período de la pandemia del COVID-19, caracterizado por restricciones a la movilidad y a las actividades productivas, que afectaron diferentes sectores de la economía, entre ellos la construcción.

Extracción señales variable 1

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

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

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

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

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

Interpretación descomposición temporal

  1. Gráfica Estacional: se observa un patrón que se repite de manera consistente cada año, con una disminución marcada en enero y los niveles más altos alrededor de julio.
  2. Gráfica Residuo: se observa que el componente residual se mantiene alrededor de cero, aunque presenta algunos valores atípicos, destacándose una caída muy pronunciada alrededor de 2020, coincidente con el período de la pandemia del COVID-19.
  3. Gráfica Tendencia: se observa un crecimiento general de los despachos durante los primeros años, seguido de algunas variaciones y una caída significativa alrededor de 2020, asociada al comportamiento extraordinario observado durante la pandemia. Posteriormente, la tendencia presenta una recuperación y alcanza niveles más elevados alrededor de 2021-2022, Luego, se observa una disminución gradual hasta 2024 y una recuperación nuevamente hacia 2025.

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

library(ggplot2)
library(plotly)

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

# Crear el gráfico
grafico_serie_var2 <- ggplot(data_col, aes(
  x = seq.Date(
    from = as.Date("2012-01-01"),
    by = "month",
    length.out = nrow(data_col)
  ),
  y = variable2
)) +
  geom_line(color = "grey", linewidth = 0.4) +
  geom_point(color = "black", size = 0.1) +
  ggtitle("V2 PNCEM: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad variable 2") +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma)

ggplotly(grafico_serie_var2)

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

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

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var2, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Descomposición temporal de la V2: PNCEM",
       x = "Tiempo",
       y = "Valor")

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

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

library(ggplot2)
library(plotly)

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

# Crear el gráfico
grafico_serie_var3 <- ggplot(data_col, aes(
  x = seq.Date(
    from = as.Date("2012-01-01"),
    by = "month",
    length.out = nrow(data_col)
  ),
  y = variable3
)) +
  geom_line(color = "grey", linewidth = 0.4) +
  geom_point(color = "black", size = 0.1) +
  ggtitle("V3 LICC: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad variable 3") +
  theme_minimal() + 
  scale_y_continuous(labels = scales::comma) 

ggplotly(grafico_serie_var3)

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

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

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Descomposición temporal de la V3: LICC",
       x = "Tiempo",
       y = "Valor")

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

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 = "Serie Original"),
    size = 0.5,
    linetype = "solid"
  ) +
  geom_line(
    aes(x = fechas_var1, y = variable1_sa, color = "Serie Ajustada"),
    size = 0.6,
    linetype = "solid"
  ) +
  scale_color_manual(
    values = c(
      "Serie Original" = "grey",
      "Serie Ajustada" = "blue"
    )
  ) +
  ggtitle("V1 - DECEM: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 1") +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)

Interpretación serie original vs serie ajustada por estacionalidad

Se observa en la gráfica anterior que ambas series mantienen un comportamiento similar en cuanto a su tendencia general, aunque la serie ajustada presenta algunas diferencias respecto a los valores originales debido a la eliminación del componente estacional. La serie ajustada conserva la caída significativa presentada alrededor de 2020, lo que indica que este comportamiento no se debe únicamente a la estacionalidad, sino que corresponde a la pandemia del Covid-19.

En términos generales, el ajuste estacional permite identificar con mayor claridad el comportamiento subyacente de los despachos, mostrando que, una vez eliminado el efecto estacional, se mantiene la tendencia general de la serie y permanecen los movimientos extraordinarios que no corresponden a patrones estacionales recurrentes..

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 = "Serie Original"),
    size = 0.5,
    linetype = "solid"
  ) +
  geom_line(
    aes(x = fechas_var2, y = variable2_sa, color = "Serie Ajustada"),
    size = 0.6,
    linetype = "solid"
  ) +
  scale_color_manual(
    values = c(
      "Serie Original" = "grey",
      "Serie Ajustada" = "red"
    )
  ) +
  ggtitle("V2 - PNCEM: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 2") +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)

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 = "Serie Original"),
    size = 0.5,
    linetype = "solid"
  ) +
  geom_line(
    aes(x = fechas_var3, y = variable3_sa, color = "Serie Ajustada"),
    size = 0.6,
    linetype = "solid"
  ) +
  scale_color_manual(
    values = c(
      "Serie Original" = "grey",
      "Serie Ajustada" = "green"
    )
  ) +
  ggtitle("V3 - LICC: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 3") +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)

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" = "blue")) +
  ggtitle("V1 - DECEM: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  theme_minimal() +
  scale_y_continuous(labels = scales::comma) +
  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)

Interpretación serie original vs tendencia

En la gráfica se observa cómo la tendencia permite suavizar las fluctuaciones mensuales presentes en la serie original y facilita la identificación del comportamiento de largo plazo de los despachos. Entre 2012 y 2016 se evidencia una trayectoria creciente, alcanzando niveles superiores. Posteriormente, entre 2017 y 2019, la tendencia presenta un comportamiento relativamente estable, con algunas variaciones.

Alrededor de 2020 se observa una disminución significativa de la tendencia, coincidente temporalmente con el período de la pandemia del COVID-19. Posteriormente, se presenta una recuperación importante entre 2021 y 2023, alcanzando niveles elevados. Durante 2024 se observa una disminución gradual y, hacia 2025 y 2026, la tendencia retoma una trayectoria creciente.

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" = "red")) +
  ggtitle("V2 - PNCEM: 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" = "green")) +
  ggtitle("V3 - LICC: 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)

Tasa de crecimiento de la serie

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**

library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
  geom_line(
    aes(
      x = fechas_corregidas_var1,
      y = tasa_crecimiento_var1,
      color = "Serie Original"
    ),
    size = 0.7
  ) +
  geom_line(
    aes(
      x = fechas_corregidas_var1,
      y = tasa_tendencia_var1,
      color = "Tendencia"
    ),
    size = 0.8,
    linetype = "dashed"
  ) +
  scale_color_manual(
    values = c(
      "Serie Original" = "grey",
      "Tendencia" = "blue"
    )
  ) +
  ggtitle("V1 - DECEM: Tasa de crecimiento anual % de la serie original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

Interpretación tasa de crecimiento anual % de la serie original y la tendencia

La gráfica muestra la variación interanual de los despachos de cemento y permite comparar el comportamiento de la serie original con el de su tendencia. Durante 2020 se observa una variación negativa significativa, coincidente con el período de la pandemia del COVID-19 y las restricciones implementadas durante la emergencia sanitaria, que afectaron la actividad económica y el sector de la construcción.

Posteriormente, en 2021 se presenta un fuerte repunte positivo en la serie original. Este comportamiento debe interpretarse teniendo en cuenta el efecto de base generado por los bajos niveles registrados durante 2020, por lo que la variación porcentual resulta excepcionalmente elevada. La tendencia también presenta una recuperación, aunque de menor magnitud debido al efecto de suavizamiento.

Después de este período, las tasas de variación de la serie original presentan fluctuaciones alrededor de cero, mientras que la tendencia muestra movimientos más moderados. Hacia 2025 y 2026 se observa una recuperación de las tasas, lo que coincide con el aumento reciente identificado en la tendencia de los despachos.

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 = "Serie Original"
    ),
    size = 0.7
  ) +
  
  geom_line(
    aes(
      x = fechas_corregidas_var2,
      y = tasa_tendencia_var2,
      color = "Tendencia"
    ),
    size = 0.8,
    linetype = "dashed"
  ) +
  
  scale_color_manual(
    values = c(
      "Serie Original" = "grey",
      "Tendencia" = "red"
    )
  ) +
  
  ggtitle("V2 - PNCEM: Crecimiento anual y tendencia (2022-2025)") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal() +
  
  # Mostrar únicamente el período posterior al pico de 2021
  coord_cartesian(
    xlim = as.Date(c("2022-01-01", "2025-12-31"))
  )


# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

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 2
grafico_crecimiento_var3 <- ggplot() +
  geom_line(
    aes(
      x = fechas_corregidas_var3,
      y = tasa_crecimiento_var3,
      color = "Serie Original"
    ),
    size = 0.7
  ) +
  geom_line(
    aes(
      x = fechas_corregidas_var3,
      y = tasa_tendencia_var3,
      color = "Tendencia"
    ),
    size = 0.8,
    linetype = "dashed"
  ) +
  scale_color_manual(
    values = c(
      "Serie Original" = "grey",
      "Tendencia" = "green"
    )
  ) +
  ggtitle("V3 - LICC: Tasa de crecimiento anual % de la serie original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()


# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

Analizar la tasa de crecimiento anual ayuda a detectar cambios en el entorno económico que afectan el sector. Se pueden prever crisis o períodos de auge y prepararse para ellos.

Modelo ARIMA

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

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

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

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

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

# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 2025 y  el conjunto de prueba o test: octubre, noviembre y diciembre de 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(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.4746  -0.9326
## s.e.  0.0876   0.0415
## 
## sigma^2 = 8.836e+09:  log likelihood = -2110.23
## AIC=4226.47   AICc=4226.62   BIC=4235.77
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 9950.195 93141.17 65425.63 -0.6702031 7.726189 0.8235653
##                     ACF1
## Training set -0.01847016

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  0.474648   0.087608   5.4179 6.031e-08 ***
## ma1 -0.932593   0.041529 -22.4564 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajuste del modelo ARIMA(1,1,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(1, 1, 1))  # Especificamos directamente (p=1, d=1, q=1)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.4746  -0.9326
## s.e.  0.0876   0.0415
## 
## sigma^2 = 8.836e+09:  log likelihood = -2110.23
## AIC=4226.47   AICc=4226.62   BIC=4235.77
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 9950.195 93141.17 65425.63 -0.6702031 7.726189 0.8235653
##                     ACF1
## Training set -0.01847016

Validación de residuales o errores del modelo

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

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 28.583, df = 22, p-value = 0.1572
## 
## Model df: 2.   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 fechas reales del período de prueba
fechas_pronostico <- seq(
  from = as.Date("2025-10-01"),
  by = "month",
  length.out = length(test_ts)
)


# 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 (1,1,1):

El modelo automático (1,1,1) presenta un comportamiento similar al observado durante el período de prueba, logrando identificar la tendencia general descendente de los despachos. Sin embargo, se evidencian diferencias entre los valores pronosticados y los observados, principalmente al inicio del período, donde el modelo presenta un mayor nivel de error. A medida que avanza el período de pronóstico, ambas series se aproximan y mantienen una trayectoria similar.

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

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

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

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

# Mostrar la tabla
print(forecast_table_auto)
##     Tiempo Observado Pronosticado
## 1 2025.750   1201024      1122785
## 2 2025.833   1112120      1094689
## 3 2025.917   1071996      1081353

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

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

# Cargar librerías necesarias
library(forecast)

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

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

# Mostrar el pronóstico completo
print(next_month_forecast_auto)
##     Tiempo Pronostico
## 1 2025.750    1122785
## 2 2025.833    1094689
## 3 2025.917    1081353
## 4 2026.000    1075024
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2026:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2026: 2026 = 1075023.58391786"

Modelo SARIMA automático

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

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

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

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

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

  • La serie presenta una tendencia no estacionaria, por lo que requiere una diferenciación no estacional.
  • El modelo incorpora dependencia temporal mediante un término autorregresivo AR(1) y un término de media móvil MA(1).
  • Existe un componente de media móvil estacional SMA(1), asociado a una periodicidad de 12 meses.
  • El modelo no requiere diferenciación estacional ni incorpora un componente autorregresivo estacional.
  • El modelo presenta un AIC de 4220,97 y un BIC de 4233,37. Estos criterios serán utilizados posteriormente para comparar su desempeño con otros modelos de pronóstico.

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(1,1,1)(0,0,1)[12] 
## 
## Coefficients:
##          ar1      ma1    sma1
##       0.4695  -0.9434  0.2123
## s.e.  0.0866   0.0401  0.0794
## 
## sigma^2 = 8.467e+09:  log likelihood = -2106.49
## AIC=4220.97   AICc=4221.22   BIC=4233.37

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(1,1,1)(0,0,1)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(1, 1, 1),  # (p,d,q) -> (1,1,1)
                seasonal = list(order = c(0, 0, 1),  # (P,D,Q) -> (0,0,1)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(1,1,1)(0,0,1)[12] 
## 
## Coefficients:
##          ar1      ma1    sma1
##       0.4695  -0.9434  0.2123
## s.e.  0.0866   0.0401  0.0794
## 
## sigma^2 = 8.467e+09:  log likelihood = -2106.49
## AIC=4220.97   AICc=4221.22   BIC=4233.37
## 
## Training set error measures:
##                   ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 9847.41 90894.46 63643.55 -0.6539708 7.538104 0.8011328
##                     ACF1
## Training set -0.03458722

Validación de residuales del modelo automatico SARIMA

En el correlograma de residuos siguiente se observa una mejora en la correlación de los residuos frente a los dos modelos anteriores. Sin embargo, al comparar los valores reales vs. pronosticados se determina una poca coincidencia, principalmente al inicio del período de prueba. Por lo anterior, aunque el modelo SARIMA presenta una mejora en el ajuste, se debe evaluar el desempeño de los pronósticos para determinar cuál modelo funciona mejor.

# 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(1,1,1)(0,0,1)[12]
## Q* = 14.332, df = 21, p-value = 0.8549
## 
## Model df: 3.   Total lags used: 24

Pronóstico con el modelo SARIMA

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

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

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

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

Interpretación modelo automatico SARIMA (1,1,1)(0,0,1):

El modelo automático (1,1,1)(0,0,1) presenta un desempeño favorable durante el período de prueba, logrando capturar la tendencia descendente de los despachos observados. Aunque inicialmente se presenta una diferencia entre los valores pronosticados y los observados, esta brecha disminuye progresivamente y hacia el final del período ambas series presentan valores muy cercanos. La incorporación del componente estacional permite considerar el comportamiento periódico de la serie, identificado previamente en el análisis de descomposición.

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

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

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

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

# Mostrar la tabla
print(forecast_table)
##     Tiempo Observado Pronosticado
## 1 2025.750   1201024      1131026
## 2 2025.833   1112120      1087910
## 3 2025.917   1071996      1071225

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

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

# Cargar librerías necesarias
library(forecast)

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

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

# Mostrar el pronóstico completo
print(next_month_forecast)
##     Tiempo Pronostico
## 1 2025.750    1131026
## 2 2025.833    1087910
## 3 2025.917    1071225
## 4 2026.000    1044283
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2026 = 1044282.74205062"

Conclusión:

El análisis de los despachos de cemento permitió identificar un comportamiento variable de la serie, con una tendencia general que presenta cambios importantes a lo largo del período analizado. Uno de los puntos más relevantes se presenta alrededor de 2020, cuando se evidencia una disminución significativa asociada al impacto de la pandemia, seguida de una recuperación posterior.

La comparación entre los modelos evaluados mostró que la incorporación del componente estacional mejora el comportamiento del pronóstico. En este sentido, el modelo SARIMA (1,1,1)(0,0,1)[12] presentó un mejor desempeño frente al modelo ARIMA, por lo que fue seleccionado como el modelo más adecuado para realizar la proyección de la variable 1.

En el período de prueba, correspondiente a octubre, noviembre y diciembre de 2025, el pronóstico logró seguir la tendencia descendente de los valores observados. Aunque existen diferencias entre ambas series, estas se reducen hacia diciembre, donde el valor pronosticado se aproxima al comportamiento observado.

A partir del modelo seleccionado, para enero de 2026 se obtiene un pronóstico de aproximadamente 1.044.283, frente a 1.071.996 registrados en diciembre de 2025, lo que representa una disminución aproximada del 2,58 %. Este resultado sugiere que, en el corto plazo, podría mantenerse una dinámica descendente en los despachos de cemento.

Desde una perspectiva empresarial, este comportamiento constituye una señal para fortalecer la planeación de producción, comercialización e inventarios, especialmente ante escenarios de menor demanda. El pronóstico puede utilizarse como una herramienta de apoyo para anticipar cambios en los despachos y ajustar oportunamente las decisiones operativas y comerciales.

En conclusión, el modelo SARIMA seleccionado permite representar de manera más adecuada el comportamiento de la serie al considerar tanto su dinámica no estacional como el componente estacional, convirtiéndose en una herramienta útil para apoyar la toma de decisiones y la planeación de los despachos de cemento.