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/Felipe/OneDrive - PUJ Cali/M Finanzas/Analitica Negocios/Modulo 2/Base Caso2.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))

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

Variable 1

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

Variable 2

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

Variable 3

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

#INTRUDUCCIÓN

El sector de la construcción en Colombia ha sido históricamente uno de los principales motores del crecimiento económico, caracterizado por un comportamiento altamente cíclico y sensible a variables macroeconómicas como las tasas de interés, el acceso al crédito, la confianza del consumidor y la política pública de vivienda.

En décadas recientes, el sector ha consolidado su importancia no solo como generador de empleo, sino también como dinamizador de industrias conexas como el cemento, el acero y los materiales de construcción. En este contexto, indicadores como los despachos de cemento (DECEM) y la producción de concreto se han convertido en variables clave para medir la actividad real del sector, mientras que las licencias de construcción (LICC) funcionan como un indicador adelantado del comportamiento futuro.

Sin embargo, el sector también ha mostrado una alta vulnerabilidad a choques externos. Un ejemplo claro fue la crisis derivada de la pandemia de COVID-19 en 2020, donde el sector registró una caída cercana al 27,7% en su PIB, acompañada de reducciones significativas en licencias de construcción (−28%), despachos de cemento (−10%) y producción de concreto (−28%), evidenciando una contracción generalizada en toda la cadena productiva.

Posteriormente, el sector mostró señales de recuperación, aunque con un comportamiento heterogéneo. En años recientes, se ha observado un panorama mixto: mientras algunos segmentos como las obras civiles han mostrado dinamismo, el subsector de edificaciones —especialmente vivienda— ha presentado desaceleraciones, reflejadas en caídas en licencias y producción de concreto.

En este contexto, el análisis conjunto de las variables LICC (licencias de construcción), DECEM (despachos de cemento) y CONCRETO (producción de concreto premezclado) permite capturar de manera integral el comportamiento del sector, desde la intención de construir (indicador adelantado), pasando por la demanda efectiva de insumos, hasta la ejecución real de obras. Este enfoque resulta particularmente relevante para identificar patrones, entender las dinámicas del sector y construir modelos de pronóstico que permitan anticipar su evolución futura y apoyar la toma de decisiones estratégicas.

Variables:

LICC (Licencias de construcción): indicador adelantado del sector (demanda futura) DECEM (Distribución de cemento): proxy de actividad constructora real CONCRETO: refleja ejecución directa de obras

Justificación: Estas tres variables capturan el ciclo completo del sector:

LICC → intención CEMENTO / CONCRETO → ejecución

¿Qué esta pasando en el sector?

Tendencia

El comportamiento del sector construcción en Colombia evidencia una dinámica cíclica estructural, con fases claras de expansión, desaceleración y choques externos, lo cual es consistente con la naturaleza del sector a nivel global.

Fase de expansión (2012–2015/2016)

Durante este periodo se observa una tendencia creciente en las licencias de construcción (LICC), acompañada por incrementos en los despachos de cemento (DECEM) y la producción de concreto, lo que refleja:

Mayor confianza del consumidor

Condiciones favorables de crédito hipotecario Políticas de estímulo a la vivienda (subsidios)

En esta fase, el crecimiento de LICENCIAS DE CONSTRUCCIÓN actúa como un indicador adelantado que impulsa posteriormente el aumento en la demanda de insumos (cemento y concreto).

Fase de desaceleración (2016–2019)

Posteriormente, el sector entra en una fase de moderación del crecimiento, evidenciada en:

Menor dinamismo en LICC Estancamiento o crecimiento más lento en DECEM Ajuste en la producción de concreto

Esto puede explicarse por:

Incremento en tasas de interés Menor crecimiento económico Reducción en la demanda de vivienda

El sector comienza a perder impulso debido a condiciones macroeconómicas menos favorables, lo que afecta principalmente la intención de nuevos proyectos.

Choque estructural (2020 – pandemia)

En este punto se presenta una ruptura abrupta en la tendencia, reflejada en:

Caídas pronunciadas en LICC Disminución significativa en DECEM Contracción fuerte en CONCRETO

Este choque afecta simultáneamente la intención, ejecución y demanda, rompiendo la dinámica normal del ciclo del sector.

Recuperación (2021–2022)

Posterior al choque, se observa una recuperación significativa, caracterizada por:

Rebote en licencias de construcción Incremento en despachos de cemento Reactivación de obras

Esto responde a:

Reactivación económica Políticas contracíclicas Demanda acumulada (efecto rebote)

Nueva desaceleración / volatilidad (2023–2025)

En los años más recientes, el sector muestra nuevamente señales de debilitamiento y alta volatilidad, evidenciadas en:

Caídas o inestabilidad en LICC Reducción en DECEM Menor dinamismo en CONCRETO

Factores explicativos:

Altas tasas de interés Inflación en costos de construcción Menor acceso a financiamiento Incertidumbre económica Políticas de Gobierno, disminución de subsidios.

El comportamiento reciente sugiere que el sector se encuentra en una fase baja del ciclo, donde la caída en licencias anticipa una posible contracción sostenida en la actividad constructora.

Estacionalidad Patrones mensuales claros: Caídas en inicio de año Recuperación en segundo semestre

Interpretación:

El sector tiene comportamiento cíclico predecible.

Componente irregular Choques visibles (ej. pandemia 2020) Volatilidad reciente

Relación entre variables

LICC lidera el ciclo CONCRETO y PNCEM reaccionan después

Interpretación:

Si caen licencias → caída futura en producción Si suben → expansión futura

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)

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

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

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)

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)

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)

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

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

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))

# Verificar longitudes
print(length(fechas_corregidas_var2))

## [1] 156

print(length(tasa_crecimiento_var2))

## [1] 156

print(length(tasa_tendencia_var2))

## [1] 156

# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("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] 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 = "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 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 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 2025 y  el conjunto de prueba o test: noviembre 2025-diciembre 2025 

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

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,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1        mean
##       0.7165  -0.4703  2044906.71
## s.e.  0.1422   0.1783    77707.53
## 
## sigma^2 = 2.947e+11:  log likelihood = -2411.45
## AIC=4830.9   AICc=4831.15   BIC=4843.33
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 2977.704 537921.3 351709.4 -6.941459 19.08833 0.6959661
##                      ACF1
## Training set -0.002851981

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        7.1651e-01  1.4215e-01  5.0405 4.643e-07 ***
## ma1       -4.7025e-01  1.7834e-01 -2.6368  0.008368 ** 
## intercept  2.0449e+06  7.7708e+04 26.3154 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo ARIMA(1,0,1) parece adecuado en el sentido de que:

todos los coeficientes son significativos; existe evidencia de dependencia temporal de corto plazo; la dinámica de la serie se explica tanto por: la influencia del valor pasado (AR1), como por la influencia de errores previos (MA1).

En lenguaje más técnico, esto sugiere que la serie tiene una estructura temporal bien capturada por un modelo mixto ARMA(1,1), ya que al tener d=0, el ARIMA(1,0,1) equivale a un ARMA(1,1).

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

# Mostrar resumen del modelo ajustado
summary(darima_auto)

## Series: train_ts 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1        mean
##       0.7165  -0.4703  2044906.71
## s.e.  0.1422   0.1783    77707.53
## 
## sigma^2 = 2.947e+11:  log likelihood = -2411.45
## AIC=4830.9   AICc=4831.15   BIC=4843.33
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 2977.704 537921.3 351709.4 -6.941459 19.08833 0.6959661
##                      ACF1
## Training set -0.002851981

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,0,1) with non-zero mean
## Q* = 34.615, df = 22, p-value = 0.04249
## 
## Model df: 2.   Total lags used: 24

La serie de residuos en el tiempo (arriba) Los residuos oscilan alrededor de cero, lo cual es positivo (cumple con media ≈ 0). Sin embargo, se observan picos muy pronunciados (outliers), especialmente alrededor de ciertos años (≈ 2015, 2020, 2023).

Interpretación:

El modelo captura la dinámica general, pero no está explicando bien choques extremos. Estos picos pueden deberse a: eventos atípicos (económicos, operativos, etc.), cambios estructurales, o necesidad de componentes adicionales (intervenciones o variables exógenas).

ACF de los residuos (abajo izquierda) La mayoría de las barras están dentro de las bandas de confianza → buena señal. Sin embargo, hay algunos rezagos significativos (por ejemplo alrededor de lag 12 y 36).

Interpretación:

Idealmente, todos los rezagos deberían estar dentro de las bandas → ruido blanco. Aquí hay evidencia de autocorrelación residual leve.

Esto sugiere que:

El modelo no captura completamente toda la dependencia temporal. Puede haber: - estructura no modelada, - o incluso un patrón estacional débil (lag 12) que no fue incluido.

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,0,1):El modelo automático (1,0,1) parece no pronosticar muy bien dentro de prueba. No hay un buen ajuste, no se capturan muy bien los puntos de quiebre. Es un modelo tentativo no adecuado para pronpostico fuera de muestra o a futuro cercano.

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   1704959      1909745
## 2 2025.833   1418504      1948062
## 3 2025.917   1937568      1975516

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    1909745
## 2 2025.833    1948062
## 3 2025.917    1975516
## 4 2026.000    1995188

# 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 = 1995187.63452649"

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

(0,1,1): Parte ARIMA no estacional: 1 términos autorregresivos (AR). 0 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,0,1)(0,0,1)[12] sugiere que:

La serie tiene una tendencia no estacionaria, corregida con cero 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 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,0,1)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    sma1        mean
##       0.6603  -0.3781  0.2685  2035603.47
## s.e.  0.1552   0.1913  0.0685    91081.87
## 
## sigma^2 = 2.716e+11:  log likelihood = -2404.65
## AIC=4819.29   AICc=4819.67   BIC=4834.82

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,0,1)(0,0,1)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(1, 0, 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,0,1)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1    sma1        mean
##       0.6603  -0.3781  0.2685  2035603.47
## s.e.  0.1552   0.1913  0.0685    91081.87
## 
## sigma^2 = 2.716e+11:  log likelihood = -2404.65
## AIC=4819.29   AICc=4819.67   BIC=4834.82
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 5711.936 514769.2 334482.9 -6.235087 18.21572 0.6618781 0.00128498

library(lmtest)

# Evaluar la significancia estadística de los coeficientes del modelo SARIMA
coeftest(auto_arima_model)

## 
## z test of coefficients:
## 
##              Estimate  Std. Error z value  Pr(>|z|)    
## ar1        6.6032e-01  1.5517e-01  4.2554 2.086e-05 ***
## ma1       -3.7814e-01  1.9133e-01 -1.9764   0.04811 *  
## sma1       2.6850e-01  6.8512e-02  3.9191 8.889e-05 ***
## intercept  2.0356e+06  9.1082e+04 22.3492 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo SARIMA (1,0,1)(0,0,1)[12] evidenció que tanto los componentes autorregresivo y de media móvil no estacionales, como el componente de media móvil estacional, fueron estadísticamente significativos. El coeficiente AR(1) positivo (φ₁ = 0.6603, p < 0.001) indica persistencia temporal en la serie, mientras que el coeficiente MA(1) negativo (θ₁ = -0.3781, p < 0.05) sugiere un efecto correctivo de los errores del periodo anterior. Adicionalmente, el componente estacional SMA(1) (Θ₁ = 0.2685, p < 0.001) confirma la presencia de una estructura cíclica con periodicidad de 12 periodos. El intercepto significativo indica que la serie fluctúa alrededor de un nivel medio estable. En conjunto, el modelo logra capturar tanto la dinámica de corto plazo como los patrones estacionales de la serie.

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 lso dos modelos anteriores. Sin embargo, al comparar los valores reales VS pronosticados se determina una poca coincidencia. Sigue funcionando mejor el modelo automatico (1,0,1)

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

El modelo SARIMA (1,0,1)(0,0,1)[12] parece mejor especificado que el ARIMA no estacional, porque:

los residuos se centran alrededor de cero; la ACF residual muestra menos evidencia de estructura remanente; la componente estacional de periodo 12 parece haber corregido parte importante de la autocorrelación.

Se evidencian dos alertas: persisten outliers importantes; la distribución residual no es perfectamente normal.

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

El modelo SARIMA logra capturar adecuadamente la tendencia y el nivel promedio de la serie, pero presenta limitaciones para anticipar cambios abruptos o choques de corto plazo, lo que se refleja en errores significativos cuando la serie experimenta caídas o recuperaciones rápidas. Esto sugiere que el modelo es adecuado para describir la dinámica estructural del sector, pero no para capturar su componente irregular o altamente volátil.

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   1704959      1763519
## 2 2025.833   1418504      1933803
## 3 2025.917   1937568      1985633

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    1763519
## 2 2025.833    1933803
## 3 2025.917    1985633
## 4 2026.000    1945390

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

## [1] "Pronóstico para enero 2026: 2026 = 1945389.91230686"

accuracy(arima_forecast_auto, test_ts)

##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set    2977.704 537921.3 351709.4  -6.941459 19.08833 0.6959661
## Test set     -257430.626 328536.5 257430.6 -17.100624 17.10062 0.5094063
##                      ACF1 Theil's U
## Training set -0.002851981        NA
## Test set     -0.592438805 0.7742387

accuracy(arima_forecast, test_ts)

##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set    5711.936 514769.2 334482.9  -6.235087 18.21572 0.6618781
## Test set     -207308.095 300706.3 207308.1 -14.080776 14.08078 0.4102234
##                     ACF1 Theil's U
## Training set  0.00128498        NA
## Test set     -0.66640872 0.7553156

Box.test(residuals(darima_auto), lag = 24, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals(darima_auto)
## X-squared = 34.615, df = 24, p-value = 0.07432

Box.test(residuals(darima), lag = 24, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals(darima)
## X-squared = 17.591, df = 24, p-value = 0.8224

library(forecast)

# Modelos
modelo_arima   <- arima(train_ts, order = c(1,0,1))
modelo_sarima  <- auto.arima(train_ts, seasonal = TRUE)

# Pronósticos
fc_arima  <- forecast(modelo_arima, h = length(test_ts))
fc_sarima <- forecast(modelo_sarima, h = length(test_ts))

acc_arima  <- accuracy(fc_arima, test_ts)
acc_sarima <- accuracy(fc_sarima, test_ts)

acc_arima

##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set    2977.704 537921.3 351709.4  -6.941459 19.08833 0.6959661
## Test set     -257430.626 328536.5 257430.6 -17.100624 17.10062 0.5094063
##                      ACF1 Theil's U
## Training set -0.002851981        NA
## Test set     -0.592438805 0.7742387

acc_sarima

##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set    5711.936 514769.2 334482.9  -6.235087 18.21572 0.6618781
## Test set     -207308.095 300706.3 207308.1 -14.080776 14.08078 0.4102234
##                     ACF1 Theil's U
## Training set  0.00128498        NA
## Test set     -0.66640872 0.7553156

comparacion <- data.frame(
  Modelo = c("ARIMA (1,0,1)", "SARIMA (1,0,1)(0,0,1)[12]"),
  RMSE = c(acc_arima[2,"RMSE"], acc_sarima[2,"RMSE"]),
  MAE  = c(acc_arima[2,"MAE"],  acc_sarima[2,"MAE"]),
  MAPE = c(acc_arima[2,"MAPE"], acc_sarima[2,"MAPE"])
)

comparacion

Modelo	RMSE	MAE	MAPE
ARIMA (1,0,1)	328536.5	257430.6	17.10062
SARIMA (1,0,1)(0,0,1)[12]	300706.3	207308.1	14.08078

Se comparó el desempeño predictivo de los modelos ARIMA (1,0,1) y SARIMA (1,0,1)(0,0,1)[12] mediante los indicadores RMSE, MAE y MAPE. Los resultados evidencian que el modelo con menor valor en estas métricas presenta mejor capacidad de pronóstico, al reducir el error entre los valores observados y estimados. No obstante, debido al tamaño reducido del conjunto de prueba, los resultados deben interpretarse con cautela.

Conclusión:

Conclusiones estratégicas 1. ¿Qué decisiones debería tomar la empresa?

La evidencia muestra que la serie presenta:

dependencia temporal (inercia del negocio) estacionalidad clara (ciclos anuales) eventos atípicos (picos abruptos)

Por tanto, la empresa debería:

✔ Planificar con enfoque estacional ✔ Implementar pronósticos como herramienta de gestión ✔ Tomar decisiones anticipadas, no reactivas El modelo permite anticipar comportamientos, no solo explicarlos. Esto reduce decisiones basadas en intuición.

¿Qué riesgos u oportunidades se identifican?

Riesgos Choques no controlados (outliers) Implica: el modelo no captura eventos extraordinarios → riesgo operativo.

Oportunidades Aprovechar la estacionalidad Los ciclos no son solo riesgo → son ventanas estratégicas

Oportunidad: campañas comerciales en meses pico optimización de precios preparación de inventarios

¿Cómo prepararse para el futuro cercano?

Adoptar un sistema de pronóstico continuo
No usar el modelo una sola vez
Actualizar mensualmente con nuevos datos

El sector está altamente expuesto a variables macroeconómicas

Los factores que explican la fase actual pueden ser a causa de:

tasas de interés altas inflación en costos menor acceso al crédito incertidumbre económica cambios en subsidios de vivienda

El sector no está cayendo por ineficiencia interna, sino por un entorno económico restrictivo.

El sector construcción en Colombia se encuentra actualmente en una fase contractiva del ciclo, caracterizada por una caída en la intención de nuevos proyectos (licencias), lo que anticipa una reducción sostenida en la actividad futura. Esta desaceleración está impulsada por condiciones macroeconómicas restrictivas, como altas tasas de interés, inflación en costos y menor acceso al financiamiento. La contracción no es aislada, sino que afecta toda la cadena productiva, desde la planificación hasta la ejecución de obras. Aunque persisten patrones estacionales, el comportamiento reciente responde principalmente a factores estructurales y no a variaciones temporales, lo que incrementa la incertidumbre y el riesgo en el corto plazo.

Los resultados confirman que la dinámica del sector presenta persistencia temporal y estacionalidad anual, por lo que un modelo SARIMA ofrece una representación más adecuada que un ARIMA no estacional. En consecuencia, la evolución reciente de las licencias, junto con la debilidad observada en los indicadores de ejecución, sugiere que el sector enfrentará en el corto plazo un entorno de menor dinamismo, mayor volatilidad y necesidad de planeación estratégica basada en pronósticos.

Series de tiempo: Extracción de señales y Modelo ARIMA - CONSTRUCTORA XYZ

Mario Felipe Montenegro

Juan Sebastián Blanco

2025-03-19

Extracción de señales

Modelo ARIMA

Modelo ARIMA automático normal (sin tener en cuenta el factor estacional)

Modelo SARIMA automático