Dinámica de la actividad constructora en Colombia: análisis de la producción de cemento, concreto premezclado y licencias de construcción

Introducción

En Colombia, el sector construcción es una de las actividades económicas mas importantes por su participación en el PIB, en generación de empleo e inversión nacional y extranjera.

Uno de los materiales indispensables en este sector es la producción de cemento, es por eso que dentro del presente articulo presentamos el comportamiento que han tenido dos variables importantes en este sector que afectan directamente la cantidad de cemento que se produce en el país, estas variables son la cantidad de metros cúbicos de concreto premezclado y el área aprobada en licencias de construcción en el país, el analisis se hace para el periodo comprendido entre enero 2012 y diciembre 2025.

Descripción y comportamiento de las variables

Para adentrar en el

****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/aleja/Downloads/Fernanda Arteaga/Analitica de datos/II Evaluación/Datos construcción.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric"))
View(data_col)

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

Variable 1

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

Variable 2

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

Variable 3

# Convertir/declarar el area de licencias de construcción aprobada en serie de tiempo mensual
variable3_ts <- ts(data_col$LICC, start = c(2012, 1), frequency = 12)

Extracción de señales

A continuación se muestre el comportamiento de las variables para el periodo de analisis que abarca desde enero 2012 a diciembre 2025.

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 = "orange", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Serie original - Producción de cemento (Toneladas)") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

library(ggplot2)
library(plotly)

# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable2_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 = "orange", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 2: Concreto premezclado (Metros cubicos)") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

library(ggplot2)
library(plotly)

# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable3_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 = "orange", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 3: Licencias de construcción (área aprobada)") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Comportamiento de las variables en conjunto:

library(ggplot2)
library(plotly)

# Convertir las 3 series a numérico
data_col$variable1_ts <- as.numeric(variable1_ts)
data_col$variable2_ts <- as.numeric(variable2_ts)
data_col$variable3_ts <- as.numeric(variable3_ts)

# Crear el gráfico
grafico_series <- ggplot(data_col, aes(
  x = seq.Date(
    from = as.Date("2012-01-01"),
    by = "month",
    length.out = nrow(data_col)
  )
)) +
  
  geom_line(aes(y = variable1_ts, color = "Ton. Producción de cemento"), linewidth = 0.6) +
  geom_line(aes(y = variable2_ts, color = "Mt. Concreto Premezclado"), linewidth = 0.6) +
  geom_line(aes(y = variable3_ts, color = "Área. Licencias de construcción"), linewidth = 0.6) +
  
  labs(
    title = "Evolución de las tres variables",
    x = "Tiempo",
    y = "Valor",
    color = "Variables"
  ) +
  theme_minimal()

ggplotly(grafico_series)

La grafica anterior permite observar que la producción de concreto premezclado y producción de cemento tienen un comportamiento similar, sin que esto signifique que la una sea causa de la otra. Es importante resaltar que aunque cada variable tiene medidas diferentes, en materia de analisis sobre el comportamiento de las mismas, el grafico presentado permite evaluar el mismo.

Sin embargo, en materia de área de licencias de construcción aprobada, el comportamiento es ligeramente mas volatil con picos bastantes altos.

Se tiene entonces y de acuerdo al grafico presentado que las tres series están expuestas a factores comunes relacionados con la dinámica económica y constructora.

Importante resaltar que para el año 2020 la caida de las 3 variables es significativa y su principal razón es la pandemia COVID 19 lo que estanco todo tipo de actividades economicas.

Extracción de señales

A continuación se hace la extracción de señales para cada una de las variables con el proposito de identificar una serie temporal en los diferentes componentes que explican su comportamiento de cada una de las variables.

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 -Producción de cemento (Ton)",
       x = "Tiempo",
       y = "Valor")

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

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
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 - Concreto premezclado (Metros cubicos)",
       x = "Tiempo",
       y = "Valor")

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

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

Tras realizar la extracción de variables se encontró que estas tienen un componente estacional con una periodicidad anual, razón por la cual es necesario realizar un ajuste estacional.

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 = "orange", 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("Producción de cemento (Ton):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)

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 = "orange", 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("Concreto premezclado (metro cubico):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 = "orange", 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("Licencias de construcción (Área aprobada):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" = "orange", "Tendencia" = "black")) +
  ggtitle("Producción de cemento: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

Tendencia Variable 2

library(ggplot2)
library(plotly)

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

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "orange", "Tendencia" = "black")) +
  ggtitle("Concreto premezclado (Metros cubicos): 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" = "orange", "Tendencia" = "black")) +
  ggtitle("Licencias de construcción (Área aprobada): 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)

Tasas de crecimiento

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

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

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

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

# Verificar longitudes
print(length(fechas_corregidas_var1))

## [1] 156

print(length(tasa_crecimiento_var1))

## [1] 156

print(length(tasa_tendencia_var1))

## [1] 156

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

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 = "orange", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Producción de cemento (Ton): 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)

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 = "orange", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Concreto premezclado (Met Cúbico): 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] 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 = "orange", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Licencias de construcción (Área aprobada): 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 diciembre de 2024. Conjunto de prueba (test): Datos desde enero 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 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) - 12 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable1_ts, end = c(2024, 12))  # Entrenamiento hasta septiembre 2024
test_ts <- window(variable1_ts, start = c(2025, 1))  # 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.3311  -0.8944
## s.e.  0.0941   0.0458
## 
## sigma^2 = 1.105e+10:  log likelihood = -2011.7
## AIC=4029.4   AICc=4029.56   BIC=4038.53
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 10757.97 104114.2 71724.91 -1.490755 8.833092 0.8967655
##                     ACF1
## Training set -0.03495811

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.331132   0.094125   3.518 0.0004348 ***
## ma1 -0.894380   0.045792 -19.532 < 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(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.3311  -0.8944
## s.e.  0.0941   0.0458
## 
## sigma^2 = 1.105e+10:  log likelihood = -2011.7
## AIC=4029.4   AICc=4029.56   BIC=4038.53
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 10757.97 104114.2 71724.91 -1.490755 8.833092 0.8967655
##                     ACF1
## Training set -0.03495811

##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* = 46.827, df = 22, p-value = 0.001546
## 
## Model df: 2.   Total lags used: 24

Box.test(
  residuals(darima_auto),
  lag = 12,
  type = "Ljung-Box",
  fitdf = 2
)

## 
##  Box-Ljung test
## 
## data:  residuals(darima_auto)
## X-squared = 29.478, df = 10, p-value = 0.001042

Los residuos del modelo ARIMA(1,1,1) todavía presentan autocorrelación significativa.

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 (1,1,1):El modelo automático (1,1,1) no es un modelo util, no es valido porque no capturan ni los puntos de quiebre y tampoco muestra el comportamiento real. Por esta razón no se avanza con el modelo.

Modelo SARIMA automático

Se procede a realizar el modela SARIMA puessto que el modelo tradicional anterior es un modelo no util. Así se espera una mejora al modelo arima tradicional ya que recoge el efecto estacional de las variables.

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.3726  -0.9185  0.3112
## s.e.  0.0947   0.0459  0.0770
## 
## sigma^2 = 9.968e+09:  log likelihood = -2003.76
## AIC=4015.52   AICc=4015.79   BIC=4027.7

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)(0,0,1)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(0, 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(0,1,1)(0,0,1)[12] 
## 
## Coefficients:
##           ma1    sma1
##       -0.6520  0.2995
## s.e.   0.1053  0.0764
## 
## sigma^2 = 1.068e+10:  log likelihood = -2009.37
## AIC=4024.74   AICc=4024.9   BIC=4033.87
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 3735.591 102345.5 64948.53 -2.061364 8.203788 0.8120415 0.1310891

##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 dos modelos anteriores*

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

Pronóstico con el modelo SARIMA

Pronóstico con el modelo SARIMA dentro del set de prueba-Gráfico líneas

# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(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  2025.000  961894.1      1105499
## 2  2025.083 1066674.8      1129406
## 3  2025.167 1200759.4      1135013
## 4  2025.250 1079403.9      1152150
## 5  2025.333 1212393.3      1128942
## 6  2025.417 1097446.2      1128285
## 7  2025.500 1262045.0      1145619
## 8  2025.583 1257260.6      1161376
## 9  2025.667 1237727.4      1130025
## 10 2025.750 1271214.1      1151660
## 11 2025.833 1226493.9      1146083
## 12 2025.917 1200229.7      1137998

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 3 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) + 3)

# 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.000    1105499
## 2  2025.083    1129406
## 3  2025.167    1135013
## 4  2025.250    1152150
## 5  2025.333    1128942
## 6  2025.417    1128285
## 7  2025.500    1145619
## 8  2025.583    1161376
## 9  2025.667    1130025
## 10 2025.750    1151660
## 11 2025.833    1146083
## 12 2025.917    1137998
## 13 2026.000    1142876
## 14 2026.083    1144693
## 15 2026.167    1145371

# 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.16666666667 = 1145370.54448311"

El modelo SARIMA (0,1,1)(0,0,1)[12] mostro mejor desempeño mostró en la comparación entre los datos reales y los pronosticados dentro del periodo de prueba (Ene, feb y marzo 2026). Destacó por su mayor precisión en la captura de los puntos de quiebre, lo que lo hace el más confiable.

No obstante, al analizar los residuos de los modelos, se identifican posibles áreas de mejora para robustecer los pronósticos en los tres casos evaluados. Algunas estrategias podrían incluir la aplicación de una transformación logarítmica o trabajar desde el inicio con la serie ajustada por estacionalidad.

modelo_sarima <- auto.arima(
  train_ts,
  seasonal = TRUE,
  stepwise = FALSE,
  approximation = FALSE
)

summary(modelo_sarima)

## Series: train_ts 
## ARIMA(1,1,1)(1,0,0)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1
##       0.3943  -0.9239  0.3348
## s.e.  0.0977   0.0484  0.0759
## 
## sigma^2 = 9.808e+09:  log likelihood = -2002.59
## AIC=4013.18   AICc=4013.45   BIC=4025.36
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 8448.329 97759.77 64864.66 -1.502661 8.054301 0.8109929
##                     ACF1
## Training set -0.03353599

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

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(1,0,0)[12]
## Q* = 11.217, df = 21, p-value = 0.9584
## 
## Model df: 3.   Total lags used: 24

# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(modelo_sarima, 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

Box.test(
  residuals(modelo_sarima),
  lag = 12,
  type = "Ljung-Box",
  fitdf = 3
)

## 
##  Box-Ljung test
## 
## data:  residuals(modelo_sarima)
## X-squared = 5.9086, df = 9, p-value = 0.749

pronostico <- forecast(modelo_sarima, h = length(test_ts))

accuracy(pronostico, test_ts)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  8448.329 97759.77 64864.66 -1.502661 8.054301 0.8109929
## Test set     38725.267 88136.10 83032.93  2.743728 7.089190 1.0381480
##                     ACF1 Theil's U
## Training set -0.03353599        NA
## Test set      0.21061642 0.8034557

El modelo SARIMA(1,1,1)(1,0,0)[12] presentó un mejor desempeño para entrenamiento y prueba. Se obtuvo un p-value de 0,749, por lo que no se encontró evidencia estadísticamente significativa de autocorrelación en los residuos.

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

# Generar pronóstico con el modelo SARIMA identificado
sarima_forecast <- forecast(modelo_sarima, h = length(test_ts))

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

# Mostrar la tabla
print(forecast_table_sarima)

##      Tiempo Observado Pronosticado
## 1  2025.000  961894.1      1086799
## 2  2025.083 1066674.8      1122037
## 3  2025.167 1200759.4      1132897
## 4  2025.250 1079403.9      1144597
## 5  2025.333 1212393.3      1127455
## 6  2025.417 1097446.2      1117832
## 7  2025.500 1262045.0      1147582
## 8  2025.583 1257260.6      1158816
## 9  2025.667 1237727.4      1132624
## 10 2025.750 1271214.1      1156921
## 11 2025.833 1226493.9      1136680
## 12 2025.917 1200229.7      1144599

# Cargar librerías necesarias
library(forecast)

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

# 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.000    1086799
## 2  2025.083    1122037
## 3  2025.167    1132897
## 4  2025.250    1144597
## 5  2025.333    1127455
## 6  2025.417    1117832
## 7  2025.500    1147582
## 8  2025.583    1158816
## 9  2025.667    1132624
## 10 2025.750    1156921
## 11 2025.833    1136680
## 12 2025.917    1144599
## 13 2026.000    1121348
## 14 2026.083    1133145
## 15 2026.167    1136780

# 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.16666666667 = 1136779.99816516"

Conclusión

El sector construcción muestra estabilidad con una suave tendencia de crecimiento. Se observa entonces que si bien para enero 2026 la proyección de producción de cemento es de una leve caida, en los meses de febrero y marzo se espera que se recupere nuevamente.

En materia de uso de modelos, inicialmente se uso el modelo ARIMA (1,1,1), sin embargo este modelo y dada la estacionalidad que presentan las variables resulto no ser viable, pues no lograba captar los puntos de quiebre y su pronostico era lineal. Por esta razón, fue necesario el uso de un model SARIMA que tuviera en cuenta la estacionalidad presentada los meses 12 de cada año quedando de la siguiente manera:

ARIMA(1,1,1)(1,0,0)[12]