Introducción

Para el caso número 2, extracción de señales y modelo ARIMA, se escoge el sector de construcción - producción de cemento y la empresa Cementos San Marcos, las variables escogidas son:

Producción de cemento: PNCEM (nacional) Despachos de cemento: DCEM (nacional) Despacho de cemento en el valle del cauca: CEM_V (Valle)

La variable principal de analisis es: la variabel 1 PNCEM

Instalar/Cargar librerias necesarias para el análisis

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

Cargar base de datos

library(readxl)
data_col <- read_excel("C:/Users/maruj/OneDrive/Escritorio/Maestria DYO/Analitica de negocio/Caso 2/Base Caso2.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric"))

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

Variable 1

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

Variable 2

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

Variable 3

# Convertir/declarar los despachos de cemento en el valle del cauca en serie de tiempo mensual
variable3_ts <- ts(data_col$CEM_V, 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("Variable 1: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Interpretación serie original

La variable de producción de cemento es una serie que no es completamente estable, se evidencia que existen subidas y bajadas frecuentes en el corto plazo.

Estas fluctuaciones pueden estar asociadas a factores estacionales (ejemplo: menor producción en algunos meses por clima o demanda o factores externos). El sector de la construcción en especial de la producción de cemento, presenta unas caidas muy fuertes en abril de 2020 producto de la pandemia, comportamiento asociado a la suspensión de la actividad productiva y obras de construcción, debido a las medidas de confinamiento y en abril de 2021 por el paro nacional, que afectaron la operación del sector en general. Para el caso en particular Cementos San Marcos, una empresa Valle Caucana, cuyas operaciones iniciaron en 2012, en sus inicios la planta tenia una capacidad de producción de 137.000 toneladas anuales, la compañia crecio rapidamente y a 2019 reporto una produccion de 650.000 toneladas año, sin embargo tambien sufrió la afectación de la pandemia y el paro nacional y para el año 2021 Cementos San Marcos reportó una producción de 600.000 toneladas Pese a las afectaciones desencadenadas por el impacto de la pandemia en el año 2020 y el paro nacional del 2021, para Cementos San Marcos fue un reto y su prioridad mantener todos sus puestos de trabajo y enfocar una estrategia que permitiera la operación de la compañía en un mercado nacional muy competido.

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(time(variable1_ts), 4),  # Tiempo repetido para cada componente (son 4 componentes)
  Value = c(stl_decomp_var1$time.series[, "seasonal"], 
            stl_decomp_var1$time.series[, "trend"], 
            stl_decomp_var1$time.series[, "remainder"], 
            variable1_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable1_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var1, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 1",
       x = "Tiempo",
       y = "Valor")+
  scale_color_manual(values = c(
    "Estacional" = "orange",
    "Tendencia" = "darkgreen",
    "Residuo" = "gray",
    "Serie Original" = "blue"
  ))


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

Interpretación descomposición temporal

En la variable de producción nacional, se evidencia en la estacionalidad, que enero es el mes donde la producción disminuye significativamente, lo cual se debe principalmente a la disminución en la actividad de la construcción, que se reduce después de la temporada de vacaciones y fin de año, así como a los factores estacionales y una posible desaceleración económica que impactan la demanda.

Por otro lado los picos fuertes de produccion son en julio de cada año y esto puede deberse a que junio es un mes de verano con la mejora de las condiciones climáticas, se intensifican los proyectos de vivienda, infraestructura y obras civiles, lo que genera una mayor demanda de cemento.

En el componente de la tendencia, muestra un crecimiento sostenido hasta 2014, de 2015 hasta finales de 2019 muestra un estancamiento sin mayores crecimientos, en 2020 nuevamente cae nuevamente producto de pandemia, el sector tiene una recuperacion importante, sin embargo en 2023 se ve afectado por la disminucion en las ventas de vivienda y las iniciaciones de nuevos proyectos, tal ves se vio afectado por el cambio de politica de otorgamiento de subsidios para la adquisicion de vivienda.

Para el residuo: La mayoría de los residuos se mantienen cerca de cero, pues la mayor parte de la variabilidad de la serie está bien explicada por tendencia + estacionalidad.

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(time(variable2_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var2$time.series[, "seasonal"], 
            stl_decomp_var2$time.series[, "trend"], 
            stl_decomp_var2$time.series[, "remainder"], 
            variable2_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable2_ts))
)

# Crear gráfico con ggplot2 con colores personalizados
p <- ggplot(stl_df_var2, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 2",
       x = "Tiempo",
       y = "Valor") +
  scale_color_manual(values = c(
    "Estacional" = "orange",
    "Tendencia" = "darkgreen",
    "Residuo" = "gray",
    "Serie Original" = "blue"
  ))

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

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(time(variable3_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_var3$time.series[, "seasonal"], 
            stl_decomp_var3$time.series[, "trend"], 
            stl_decomp_var3$time.series[, "remainder"], 
            variable3_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable3_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 3",
       x = "Tiempo",
       y = "Valor")+
  scale_color_manual(values = c(
    "Estacional" = "orange",
    "Tendencia" = "darkgreen",
    "Residuo" = "gray",
    "Serie Original" = "blue"
  ))

# 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 = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 1:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)

Interpretación serie original vs serie ajustada por estacionalidad

Podemos observar que, tras retirar la estacionalidad, la dinámica de la producción no cambia de manera estructural, pero se eliminan las oscilaciones regulares que podían ocultar la tendencia.

La caída drástica alrededor del 2020 y 2021 se mantiene en la serie ajustada, lo que nos podria confirmar que estos eventos no son estacionales, sino un choque real en la producción (pandemia y paro nacional)

Gráfico serie original VS ajustada Variable 2

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

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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)

Gráfico serie original VS ajustada Variable 3

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

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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)

Ahora graficamos serie original vs tendencia

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

Tendencia Variable 1

library(ggplot2)
library(plotly)

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

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
  geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 1: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

Interpretación serie original vs tendencia

La serie muestra que al inicio hubo crecimiento, luego se estabilizó y en 2020 tuvo una caída muy fuerte (seguramente por la pandemia). Después se recuperó, pero desde 2023 la tendencia va bajando, lo que refleja una desaceleración en la producción, principalmente por el cambio en la politica de otorgamientos de subsidios de vivienda por parte del estado y la estabilizacion de las tasas de creditos bancarios, sin embargo no podemos concluir si la tendencia continue a la baja.

Tendencia Variable 2

library(ggplot2)
library(plotly)

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

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

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

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

Tendencia Variable 3

library(ggplot2)
library(plotly)

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

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

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

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

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] 144

print(length(tasa_crecimiento_var1))

## [1] 144

print(length(tasa_tendencia_var1))

## [1] 144

*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 = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable1: 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

Bajo la tasa de crecimiento anual podemos observar que la tendencia muestra un crecimiento moderado y relativamente estable hasta 2019, la produccion de cemento venia con un rito sostenido.

Hay una caida muy significativa en 2020, producto de la pandemia y debido al confinamiento se detuvieron muchas obras; posterior en el 2021 observamos una reactivación, generando un crecimiento bastante significativo, sin embargo la tendencia empieza a descender , lo que nos indica que la producción de cemento se esta ajustando despues del efecto rebote y pareciera que vuelve a la producción normal que tenia antes de pandemia.

Para 2024 y 2025 la línea esta muy cerca del 0%, se podria inferir que el crecimiento se ha frenado. No hay una contracción fuerte, pero tampoco hay señales de expansión significativa.

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] 144

print(length(tasa_crecimiento_var2))

## [1] 144

print(length(tasa_tendencia_var2))

## [1] 144

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

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] 144

print(length(tasa_crecimiento_var3))

## [1] 144

print(length(tasa_tendencia_var3))

## [1] 144

# Gráfico de la tasa de crecimiento anual variable 3
grafico_crecimiento_var3 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable3: 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 2024. Conjunto de prueba (test): Datos desde octubre de 2024 hasta diciembre de 2024.

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

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

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

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

Modelo ARIMA automático normal (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.3336  -0.8940
## s.e.  0.0957   0.0472
## 
## sigma^2 = 1.13e+10:  log likelihood = -1974.44
## AIC=3954.88   AICc=3955.04   BIC=3963.95
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 10449.67 105259.9 72659.69 -1.566531 8.971488 0.9023438
##                     ACF1
## Training set -0.03433432

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.333617   0.095680   3.4868 0.0004888 ***
## ma1 -0.893984   0.047222 -18.9315 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Ajuste del modelo ARIMA(4,1,2) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts, 
                order = c(2, 1, 1))  # Especificamos directamente (p=4, d=1, q=2)  #ojo cambiar el orden

# Mostrar resumen del modelo ajustado
summary(darima_auto)

## Series: train_ts 
## ARIMA(2,1,1) 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.3301  0.0743  -0.9112
## s.e.  0.0937  0.0894   0.0469
## 
## sigma^2 = 1.132e+10:  log likelihood = -1974.09
## AIC=3956.19   AICc=3956.46   BIC=3968.28
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 11261.48 105017.8 72027.56 -1.486761 8.910511 0.8944935
##                     ACF1
## Training set -0.01367176

Validación de residuales o errores del modelo

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

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)
## Q* = 42.087, df = 21, p-value = 0.004103
## 
## Model df: 3.   Total lags used: 24

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

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

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

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

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

Interpretación modelo automatico (2,1,1)

El modelo automático (2,1,1) parece no pronosticar dentro de prueba. No se capturan los puntos de quiebre. No es modelo adecuado para pronpostico fuera de muestra o a futuro.

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

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

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

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

# Mostrar la tabla
print(forecast_table_auto)

##     Tiempo Observado Pronosticado
## 1 2024.750   1180815      1133340
## 2 2024.833   1121219      1141154
## 3 2024.917   1145357      1143761

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

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

# Cargar librerías necesarias
library(forecast)

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

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

# Mostrar el pronóstico completo
print(next_month_forecast_auto)

##     Tiempo Pronostico
## 1 2024.750    1133340
## 2 2024.833    1141154
## 3 2024.917    1143761
## 4 2025.000    1144631

# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))

## [1] "Pronóstico para enero 2025: 2025 = 1144630.62871857"

Modelo SARIMA

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(0,1,1)(0,0,1)[12] sugiere que:

La serie tiene una tendencia no estacionaria, corregida con una diferenciación.
Existe una influencia significativa del error pasado (MA(1)).
Hay un componente estacional autorregresivo fuerte cada 12 períodos.
El ajuste es adecuado según los criterios AIC y BIC, pero se podría comparar con otros modelos para mejorar la predicción.

Identificación automática del modelo SARIMA

# Identificación automática modelo SARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)

## Series: train_ts 
## ARIMA(1,1,1)(0,0,1)[12] 
## 
## Coefficients:
##          ar1      ma1    sma1
##       0.3748  -0.9184  0.3120
## s.e.  0.0966   0.0477  0.0782
## 
## sigma^2 = 1.02e+10:  log likelihood = -1966.72
## AIC=3941.45   AICc=3941.72   BIC=3953.54

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(1, 1, 1),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(0, 0, 1),  # (P,D,Q) -> (1,0,0)
                                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.3748  -0.9184  0.3120
## s.e.  0.0966   0.0477  0.0782
## 
## sigma^2 = 1.02e+10:  log likelihood = -1966.72
## AIC=3941.45   AICc=3941.72   BIC=3953.54
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE    MAPE      MASE        ACF1
## Training set 9371.476 99667.25 67628.01 -1.49445 8.36446 0.8398565 -0.03621159

Validación de residuales del modelo automatico SARIMA

En el correlograma de residuos siguiente se observa que, mejora la correlación de los residuos frente a los modelos anteriores. Al comparar los valores reales VS pronosticados se determina que tiene coincidencia, parece pronosticar mejor dentro de prueba; hay un sobre ajuste, pero se capturan muy bien los puntos de quiebre. Es un modelo tentativo adecuado para pronpostico fuera de muestra o a futuro, ´por lo tanto este modelo funciona mejor que el modelo ARIMA.

# 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* = 15.553, df = 21, p-value = 0.7942
## 
## Model df: 3.   Total lags used: 24

Análisis del modelo SARIMA

En el autocorrelograma (ACF), las barras están dentro de las bandas de confianza (líneas punteadas azules), lo que indica que no hay autocorrelación significativa.

Pronóstico con el modelo SARIMA

Dentro de prueba

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

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

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

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

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

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

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

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

# Mostrar la tabla
print(forecast_table)

##     Tiempo Observado Pronosticado
## 1 2024.750   1180815      1147233
## 2 2024.833   1121219      1117421
## 3 2024.917   1145357      1171977

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

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 (2 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)

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

# Mostrar el pronóstico completo
print(next_month_forecast)

##     Tiempo Pronostico
## 1 2024.750    1147233
## 2 2024.833    1117421
## 3 2024.917    1171977
## 4 2025.000    1101935

# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))

## [1] "Pronóstico para enero 2025: 2025 = 1101934.82625893"

Conclusión

El modelo automático SARIMA (1,1,1)(0,0,1) fue el que mejor desempeño mostró en la comparación entre los datos reales y los pronosticados dentro del periodo de prueba (oct.nov.dic2024), se destaca su precisión en la captura de los puntos de quiebre, lo que lo hace el más confiable.

El pronóstico realizado con el modelo SARIMA, en línea con el gráfico de la tasa de crecimiento anual de la variable 1, evidencia una tendencia decreciente en la producción de cemento. Esta caída puede explicarse principalmente por la desaceleración que atraviesa el sector de la construcción en Colombia.

Dic-2024 real 1145356,506265 vs pronostico futuro ene-2025 1101934.8262, decrecimiento -3.79%

Factores como las altas tasas de interés, que encarecen los créditos hipotecarios y empresariales, la reducción de subsidios de vivienda y una inflación persistente que incrementa los costos de materiales e insumos.

Todo esto hace que la demanda de cemento baje, porque hay menos proyectos y menos compra de vivienda. Al mismo tiempo, en el lado de la oferta, producir cemento se ha vuelto más costoso por el aumento de los precios de insumos y transporte

Para empresas como Cementos San Marcos, este escenario representa un reto significativo, la menor demanda de proyectos de infraestructura y vivienda se traduce directamente en una reducción en la venta de cemento. Además, el aumento de los costos operativos presiona los márgenes de rentabilidad.

Lo anterior limita la capacidad de crecimiento de Cementos San Marcos y lo obligan a replantear estrategias comerciales y de eficiencia para sostener su posición en un mercado cada vez más competitivo y exigente.

Por lo anterior se recomienda:

Diversificar clientes: no depender solo de clientes como las constructores de viviendas, si no tambien enfocarse en obras de infrastructura vial.
Buscar plan de expansión en zonas donde aun haya dinamismo de crecimiento de construcción.
Optimización de procesos productivos, como por ejemplo el uso de soluciones energeticas o mejoras tecnológicas para reducir consumo de energía.
Uso de materias primas locales: para disminuir dependencia de insumos importados más costosos.
Mejor logística y distribución: reducir costos de transporte aprovechando su ubicación estratégica (Valle del Cauca).

Sector Construcción - Producción de Cemento

Deisy Ortega

2025-08-23