Introducción

CEMEX COLOMBIA S.A.

El sector de materiales de construcción en Colombia ha sido un motor clave para la economía del país en los últimos años, con un impacto significativo en el empleo, la inversión y el desarrollo de infraestructura. Sin embargo, ha enfrentado diversos desafíos y fluctuaciones debido a factores económicos, políticos y sociales.

CEMEX Colombia S.A., es una de las principales compañías del sector en el país. Como parte de CEMEX, una multinacional mexicana con presencia en más de 50 países, la empresa se especializa en la producción, distribución y comercialización de cemento, concreto premezclado y agregados.

Desde su llegada a Colombia, CEMEX ha sido un actor clave en el desarrollo de infraestructura y construcción de vivienda, contribuyendo al crecimiento económico y social del país. Su compromiso con la innovación y la sostenibilidad la ha llevado a implementar procesos eficientes y tecnologías avanzadas para minimizar su impacto ambiental.

Con plantas de producción en varias regiones del país y una amplia red de distribución, CEMEX Colombia atiende las necesidades de diversos sectores, incluyendo la construcción de viviendas, carreteras, puentes y grandes proyectos de infraestructura. Además, la empresa mantiene un fuerte enfoque en la responsabilidad social y el desarrollo comunitario, promoviendo iniciativas de educación, sostenibilidad y bienestar en las comunidades donde opera.

Gracias a su experiencia, calidad y compromiso con la sostenibilidad, CEMEX Colombia se ha consolidado como una empresa líder en la industria de materiales de construcción en el país.

Cargar base de datos para correlación de variables

library(readxl)
data_col <- read_excel("~/MAESTRIA FINANZAS/ANALITICA DE DATOS/Base Caso2.xlsx", col_types = c("date", 
    "numeric", "numeric", "numeric"))

1) Extracción de señales o descomposición de series temporales

En el análisis se evaluaran las siguientes variables:

-Producción de cemento - toneladas -Despacho de cemento - toneladas -Producción de concreto premezclado - metros cúbicos

Estas tres variables son fundamentales para el desarrollo económico y de infraestructura en Colombia, ya que impactan directamente sectores clave como la construcción de viviendas, obras civiles y proyectos de infraestructura vial.

Correlación entre las variables seleccionadas

Correlación entre producción de cemento y despacho de cemento

cor(data_col$PNCEM, data_col$DECEM)
## [1] 0.9342798

Correlación positiva fuerte, lo que indica que cuando la producción de cemento aumenta, los despachos también crecen casi en la misma proporción.

Correlación entre despachos de cemento y producción de concreto premezclado

cor(data_col$DECEM, data_col$CONCRETO)
## [1] 0.5936762

El coeficiente de correlación de 0.59 indica una relación positiva moderada-alta entre los despachos de cemento y la producción de concreto premezclado, es decir que a medida que los despachos de cemento aumenta, también lo hace la producción de concreto premezclado, pero no de manera totalmente proporcional.

Correlación entre producción de cemento y producción de concreto premezclado

cor(data_col$PNCEM, data_col$CONCRETO)
## [1] 0.5195142

El coeficiente de correlación de 0.51 indica una relación positiva moderada entre la producción de cemento y la producción de concreto premezclado.A medida que la producción de cemento aumenta, también lo hace la producción de concreto premezclado, pero no de manera completamente proporcional.

En el estudio de series temporales se identifican patrones subyacentes en los datos, como tendencias, estacionalidad y componentes irregulares, lo que facilita una mejor comprensión del comportamiento de las variables a analizar. Para ello, se emplearán técnicas estadísticas o métodos de descomposición

2) Pronóstico de la variable elegida para estudio

-Producción de cemento - toneladas

Para el pronóstico de esta variable se empleará el Modelo ARIMA

A lo largo del informe, se presentarán los resultados obtenidos y se discutirán sus implicaciones en el contexto empresarial de Colombia, con el fin de extraer información clave que pueda contribuir a la toma de decisiones estratégicas.

El periodo de análisis es desde el 01 de enero de 2012 hasta 01 de diciembre 2024.

Análisis exploratorio de las variables

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

Cargar base de datos

library(readxl)
data_col <- read_excel("~/MAESTRIA FINANZAS/ANALITICA DE DATOS/Base Caso2.xlsx", col_types = c("date", 
    "numeric", "numeric", "numeric"))

Un paso indispensable es declarar la variable PNCEM como serie temporal:

# Convertir/declarar el número de toneladas en serie de tiempo mensual
PNCEM_ts <- ts(data_col$PNCEM, start = c(2012, 1), frequency = 12)
PNCEM_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012  868474.0  865408.0  998847.0  852138.0  919675.0  906243.0  910743.0
## 2013  803846.0  823731.0  930060.0  888046.0  929099.0  879965.0  966824.0
## 2014  821445.8  908543.1 1089339.5 1028889.7 1074107.8  986800.1 1081296.1
## 2015  899518.5 1034033.4 1098526.0 1056733.9 1128654.1 1023576.6 1108864.4
## 2016  968000.7 1076538.1 1025163.3 1108543.4 1054354.4 1012287.1  912684.1
## 2017  897063.5 1025705.0 1099919.7  976913.7 1021787.4 1006870.5 1079282.2
## 2018  906920.5  955449.3 1066143.4 1034777.6 1007110.9  995956.4 1021557.8
## 2019  928828.6  993333.3 1147660.3  971013.0 1119226.6 1025805.1 1118126.7
## 2020 1045240.9 1042244.0  850163.8  198924.9  778730.3  977660.3 1131331.2
## 2021 1034344.7 1139862.9 1200339.1 1135881.8  880396.4 1108427.0 1201721.7
## 2022 1040160.8 1133914.0 1341585.4 1214953.5 1251173.5 1189314.5 1200980.9
## 2023 1047479.1 1140751.0 1231072.6 1111785.8 1237255.1 1116485.4 1221109.3
## 2024 1006077.9 1091777.0 1114426.0 1146361.4 1093702.8 1069209.8 1155668.7
##            Aug       Sep       Oct       Nov       Dec
## 2012  942864.0  894974.0  885860.0  931125.0  948435.0
## 2013  941125.0 1002957.0 1048557.0 1020409.0 1017314.0
## 2014 1068933.6 1071902.7 1074163.9 1069628.5 1127320.3
## 2015 1139186.2 1125937.9 1176890.0 1045135.7 1209660.3
## 2016 1130482.4 1023568.4 1076104.9 1009009.2 1098666.2
## 2017 1034803.5 1025617.1 1081553.1 1028965.7 1020361.5
## 2018 1097724.3 1083427.0 1107128.4 1094297.7 1089381.0
## 2019 1175004.6 1093948.0 1135813.5 1112915.0 1191117.6
## 2020 1095590.7 1148007.1 1190414.4 1213983.9 1176276.7
## 2021 1163338.7 1229149.7 1234917.7 1201226.8 1267497.4
## 2022 1233964.5 1279462.0 1249315.4 1193547.5 1283041.9
## 2023 1178826.2 1261109.0 1261911.4 1112027.2 1291078.1
## 2024 1189179.0 1109916.7 1180815.2 1121219.4 1145356.5

Un paso indispensable es declarar la variable DECEM como serie temporal:

# Convertir/declarar el número de toneladas en serie de tiempo mensual
DECEM_ts <- ts(data_col$DECEM, start = c(2012, 1), frequency = 12)
DECEM_ts
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012  823283.7  846615.1  950452.9  789542.0  904690.6  879219.3  879632.6
## 2013  797141.6  826921.1  808108.3  917850.5  901645.6  864810.7  980476.6
## 2014  810402.5  928329.4 1041472.4  956972.4 1034413.3  918642.5 1059514.7
## 2015  927848.3  986176.4 1079949.1  999411.6 1061957.2 1012089.6 1152372.0
## 2016  939138.8 1045943.9 1007468.5 1058103.1 1005411.0  997177.5  924987.7
## 2017  913190.4 1007968.3 1083268.2  900004.9  999046.0  983520.3 1041467.8
## 2018  909395.5  960500.3  978879.9 1030819.5  984012.0  943219.3  990012.1
## 2019  917158.5  972510.7 1035824.3  991119.7 1054963.4  979305.2 1113127.7
## 2020  994701.3 1024532.8  739689.4  242413.7  705921.0  904954.8 1093059.4
## 2021  989097.3 1079241.9 1189013.8 1050823.9  815891.2 1071880.6 1131975.9
## 2022  963889.4 1106900.3 1257124.9 1089208.1 1103450.4 1117050.3 1124719.6
## 2023  930371.1 1057742.2 1162546.0  987450.8 1126139.4 1033032.1 1051542.1
## 2024  896617.8 1032284.3  939625.7 1074893.5 1003338.7  942870.1 1050528.2
##            Aug       Sep       Oct       Nov       Dec
## 2012  906463.0  862746.7  906790.6  915069.6  831485.2
## 2013  870423.1  981342.5 1032773.3  978662.3  905880.8
## 2014 1020250.3 1085842.5 1089667.1 1040586.8  984139.6
## 2015 1108117.5 1142760.7 1161600.8 1067607.3 1106888.9
## 2016 1089617.9 1015052.3  993667.3 1017490.6 1006878.6
## 2017 1033008.5 1022836.0 1029658.7 1021276.7  948589.1
## 2018 1061280.2 1030848.1 1083220.2 1064519.7  973383.7
## 2019 1099300.5 1078411.7 1110885.7 1083302.8 1079406.9
## 2020 1052074.0 1131203.4 1169113.2 1121991.0 1054316.6
## 2021 1097788.3 1163918.5 1150044.1 1166822.2 1134856.2
## 2022 1190049.8 1189530.8 1123078.6 1118150.0 1122246.0
## 2023 1103019.7 1131366.0 1113741.8 1103166.9 1093586.7
## 2024 1074200.5 1003766.7 1073791.9 1015636.2 1002034.8

Un paso indispensable es declarar la variable CONCRETO como serie temporal:

# Convertir/declarar el número de toneladas en serie de tiempo mensual
CONCRETO_ts <- ts(data_col$CONCRETO, start = c(2012, 1), frequency = 12)
CONCRETO_ts
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2012 526135.5 584896.7 634905.0 550290.4 639649.5 620337.5 621040.8 640026.3
## 2013 553389.8 601934.0 577620.9 663616.0 668222.6 646338.8 709693.1 664980.0
## 2014 581032.9 663333.4 709935.8 692234.9 752105.4 649168.4 786839.1 698542.2
## 2015 588921.1 689289.9 735713.6 709793.6 771679.8 741080.6 834848.6 772132.9
## 2016 596128.0 733812.5 685446.2 738887.0 689904.8 721605.8 628366.5 697717.5
## 2017 537905.3 630653.8 695564.2 569457.8 629803.1 607064.4 594304.3 621530.5
## 2018 500863.0 567671.1 574867.9 589132.2 589787.6 590334.9 588338.9 621400.5
## 2019 518830.1 621271.6 655183.7 581143.7 643375.5 578861.6 651512.4 631085.4
## 2020 474608.7 575767.6 404111.7  30559.2 298568.0 439428.8 513111.8 509440.8
## 2021 391223.4 505824.3 572461.8 457335.7 398756.9 485677.2 544203.4 547982.0
## 2022 560312.2 680482.4 741457.3 653098.7 684290.7 672862.1 690243.0 760621.3
## 2023 569093.4 686987.3 743434.9 626058.5 731515.5 696036.4 683440.9 722865.7
## 2024 536390.1 676447.7 628717.1 680562.9 673810.8 615566.0 679577.3 681989.0
##           Sep      Oct      Nov      Dec
## 2012 618164.0 636544.0 613992.7 540217.9
## 2013 688418.0 737172.5 664587.9 610667.7
## 2014 743102.9 740833.7 662678.4 650640.8
## 2015 810193.7 787664.6 725788.4 705681.8
## 2016 687122.5 663349.2 635103.6 626942.7
## 2017 609339.2 616166.9 612336.1 553678.9
## 2018 610543.8 641974.6 611400.8 551577.8
## 2019 605898.1 624728.9 585884.5 548756.3
## 2020 534261.6 539360.1 483381.0 457110.3
## 2021 588808.4 569810.6 563315.5 545342.7
## 2022 752821.1 728060.0 738940.1 707833.2
## 2023 749036.2 700717.6 715122.4 641701.1
## 2024 632938.1 682140.0 638772.9 608304.9
# Calcular estadísticas descriptivas básicas
descriptive_stats <- data.frame(

  Media = mean(PNCEM_ts),
  Mediana = median(PNCEM_ts),
  DesviacionEstandar = sd(PNCEM_ts),
  CoefVar = sd(PNCEM_ts) / mean(PNCEM_ts)
)
print(descriptive_stats)
##     Media Mediana DesviacionEstandar   CoefVar
## 1 1062721 1076321           133962.9 0.1260566

Ahora si calculamos estadísticas descriptivas

En la siguiente tabla se puede observar la variable PNCEM:

  • En promedio, se ha producido alrededor de 1.062.721 de toneladas. Este es un valor de referencia para entender el comportamiento típico.

  • La mediana indica que la mitad de toneladas producidas en el periodo son de 1.076.321 y como esta cercana a la media, lo que refiere que la distribución no está muy sesgada.

  • Existe una variabilidad estandar en la cantidad de cemento producido por mes. Un valor de 117.009 indica que los datos no estan muy dispersos.

  • la variabilidad relativa Un 12.61% implica que, aunque hay oscilaciones, los datos no son extremadamente dispersos respecto al promedio, lo que sugiere una distribución relativamente estable.

# Calcular estadísticas descriptivas básicas
descriptive_stats <- data.frame(

  Media = mean(DECEM_ts),
  Mediana = median(DECEM_ts),
  DesviacionEstandar = sd(DECEM_ts),
  CoefVar = sd(DECEM_ts) / mean(DECEM_ts)
)
print(descriptive_stats)
##     Media Mediana DesviacionEstandar   CoefVar
## 1 1009828 1023684           117009.1 0.1158702

En la siguiente tabla se puede observar la variable DECEM:

  • En promedio, se ha despachado alrededor de 1.009.828 toneladas de cemento. Este es un valor de referencia para entender el comportamiento típico.

  • La mediana indica que la mitad de despachos realizados en el periodo son de 1.023.684 y como esta cercana a la media, lo que refiere que la distribución no está muy sesgada.

  • Existe una variabilidad estandar en la cantidad de despachos de cemento por mes. Un valor de 117.009 indica que los datos no fluctúan bastante alrededor del promedio.

  • la variabilidad relativa Un 11.59% implica que, aunque hay oscilaciones, los datos no son extremadamente dispersos respecto al promedio, lo que sugiere una distribución relativamente estable.

# Calcular estadísticas descriptivas básicas
descriptive_stats <- data.frame(

  Media = mean(CONCRETO_ts),
  Mediana = median(CONCRETO_ts),
  DesviacionEstandar = sd(CONCRETO_ts),
  CoefVar = sd(CONCRETO_ts) / mean(CONCRETO_ts)
)
print(descriptive_stats)
##      Media  Mediana DesviacionEstandar   CoefVar
## 1 626188.4 630869.6           100036.7 0.1597549

En la siguiente tabla se puede observar la variable CONCRETO:

  • En promedio, se ha producido alrededor de 626.188 de concreto en metros cúbicos. Este es un valor de referencia para entender el comportamiento típico.

  • La mediana indica que la mitad de metros cúbicos producidos en el periodo son de 630.870 y como esta cercana a la media, lo que refiere que la distribución no está muy sesgada.

  • Existe una variabilidad estándar en la cantidad de concreto producido por mes. Un valor de 100.036 indica que los datos no estan muy dispersos.

  • la variabilidad relativa Un 15.97% implica que, aunque hay oscilaciones, los datos no son extremadamente dispersos respecto al promedio, lo que sugiere una distribución relativamente estable.

Extracción de señales

En este análisis,la descomposición de series de tiempo permitirá comprender mejor los datos, mejorar predicciones y tomar decisiones más estratégicas. Es una herramienta clave en la analítica de negocios, especialmente en entornos donde las fluctuaciones en los datos pueden afectar inversiones, políticas económicas y estrategias empresariales.

Gráfico 1. Descomposición de la serie temporal PNCEM

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

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

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_PNCEM <- data.frame(
  Time = rep(time(PNCEM_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_PNCEM$time.series[, "seasonal"], 
            stl_decomp_PNCEM$time.series[, "trend"], 
            stl_decomp_PNCEM$time.series[, "remainder"], 
            PNCEM_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(PNCEM_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_PNCEM, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Gráfico 1. Descomposición del número de toneladas de producción de cemento",
       x = "Tiempo",
       y = "Valor")

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

En el gráfico se observa un fuerte efecto estacional en la producción de cemento, se observan picos altos y bajos en los meses de septiembre y enero respectivamente de cada año. Lo que indica que en septiembre la demanda de cemento es mayor y en enero se presenta una actividad menor de la industria lo que conlleva a una menor producción.

En la serie original se observan fluctuaciones entre los años 2012 y 2024, En el año 2020 se presenta un periodo de contracción. La serie original contiene variaciones estacionales y aleatorias que pueden dificultar la identificación de patrones reales. Por tanto, la tendencia sugiere menor volatilidad y se identifica una recuperación posterior en el año 2021 y una leve caída en el año 2024.

Adicionalmente, a través del análisis de los residuos de esta variable se puede observar un evento o situación extraordinaria en el año 2020 debido a la pandemia por COVID-19 lo que afectó la producción de cemento.

Gráfico 2. Descomposición de la serie temporal DECEM

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

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

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_DECEM <- data.frame(
  Time = rep(time(DECEM_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_DECEM$time.series[, "seasonal"], 
            stl_decomp_DECEM$time.series[, "trend"], 
            stl_decomp_DECEM$time.series[, "remainder"], 
            DECEM_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(DECEM_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_DECEM, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Gráfico 2. Descomposición del número de despachos de cemento",
       x = "Tiempo",
       y = "Valor")

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

Como se puede observar en el gráfico 2, la estacionalidad es fuerte en los despachos de cemento y tienen un patrón recurrente predecible; lo que sugiere que ciertos meses o estaciones son más activas en la construcción. Las caídas e incrementos del componente estacional ocurren en intervalos similares en los meses de enero y de julio respectivamente de cada año, lo que indica que este componente estacional se mantiene a lo largo del tiempo.
Para el componente de la tendencia se puede evidenciar un crecimiento gradual hasta 2019, seguido de una fuerte caída en 2020, tras la recuperación en 2021-2022, la tendencia muestra signos de disminución en 2024.

En el caso del residuo hay factores externos afectando los despachos de cemento, en el año 2020 se ve reflejado el impacto de la pandemia y para el año 2021 el paro nacional.

Gráfico 3. Descomposición de la serie temporal CONCRETO

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

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

# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_CONCRETO <- data.frame(
  Time = rep(time(CONCRETO_ts), 4),  # Tiempo repetido para cada componente
  Value = c(stl_decomp_CONCRETO$time.series[, "seasonal"], 
            stl_decomp_CONCRETO$time.series[, "trend"], 
            stl_decomp_CONCRETO$time.series[, "remainder"], 
            CONCRETO_ts),
  Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(CONCRETO_ts))
)

# Crear gráfico con ggplot2
p <- ggplot(stl_df_CONCRETO, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Gráfico 3. Descomposición del número de métros cúbicos de producción de concreto premezclado",
       x = "Tiempo",
       y = "Valor")

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

Como se puede identificar en el gráfico 3, la estacionalidad es fuerte en la producción de concreto premezclado y al igual que los despachos de cemento tienen un patrón recurrente predecible;Las caídas y los incrementos también ocurren en los meses de enero y de julio de cada año,lo cual es un factor clave dado que el cemento es un insumo necesario para la producción de concreto y si hay incrementos o caídas en los despachos de cemento se puede anticipar que hay un crecimiento o decrecimiento en la producción del concreto.

En la tendencia indica una caída importante en 2020, una recuperación en 2021-2022, pero con signos recientes de desaceleración.

En general los residuos son relativamente pequeños, pero hay un pico negativo pronunciado en 2020, lo que indica un evento atípico (pandemia) lo que afectó significativamente la producción.

Se crean las variables ajustada por estacionalidad

# Extraer los componentes de la descomposición
PNCEM_sa <- PNCEM_ts - stl_decomp_PNCEM$time.series[, "seasonal"]
# Extraer los componentes de la descomposición
DECEM_sa <- DECEM_ts - stl_decomp_DECEM$time.series[, "seasonal"]
# Extraer los componentes de la descomposición
CONCRETO_sa <- CONCRETO_ts - stl_decomp_CONCRETO$time.series[, "seasonal"]

Gráfico serie original vs serie ajustada variable PNCEM

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

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_PNCEM <- ggplot() +
  geom_line(aes(x = fechas_PNCEM, y = PNCEM_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_PNCEM, y = PNCEM_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Gráfico 4. Producción de cemento:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("número de toneladas") +
  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_PNCEM)

En el gráfico 4 la serie ajustada suaviza los efectos estacionales, facilitando el análisis de la tendencia subyacente. En general se observa una tendencia creciente hasta 2019, luego una fuerte caída en 2020 y una posterior recuperación con cierta estabilización en los últimos años.

La producción de cemento ha crecido en el tiempo, pero con fluctuaciones y la crisis del 2020 tuvo un impacto fuerte y claro en la producción de cemento ocasionando una ligera desaceleración.

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

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_DECEM <- ggplot() +
  geom_line(aes(x = fechas_DECEM, y = DECEM_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_DECEM, y = DECEM_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Gráfico 5. Despachos de cemento:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Número de despachos") +
  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_DECEM)

En el gráfico 5 se observa que la serie original presenta un patrón cíclico, típico de series con estacionalidad.

con relación a la serie ajustada se puede evidenciar que antes del 2020 los despachos de cemento aumentaban en el tiempo. Sin embargo, durante la pandemia se presentó una caída abrupta en la producción, seguida de una recuperación y posteriormente una desaceleración en los últimos años.

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

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_CONCRETO <- ggplot() +
  geom_line(aes(x = fechas_CONCRETO, y = CONCRETO_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_CONCRETO, y = CONCRETO_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Gráfico 6. Producción de concreto premezclado:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("metros cúbicos") +
  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_CONCRETO)

En el gráfico 6 en la serie ajustada Se observa una tendencia creciente hasta 2018, seguida de cierta estabilización. Adicionalmente, la recuperación después de la caída de 2020 parece haber llevado la producción de concreto a un nivel similar al de los años cercanos 2015 - 2016, sin una clara tendencia alcista reciente.

Extraer la tendencia PNCEM

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
PNCEM_vec <- as.numeric(PNCEM_ts)
tendencia_PNCEM <- as.numeric(stl_decomp_PNCEM$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_PNCEM <- ggplot() +
  geom_line(aes(x = fechas, y = PNCEM_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_PNCEM, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle(" Gráfico 7 PNCEM: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  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_PNCEM)

En la gráfica 7 la tendencia muestra un crecimiento en los despachos de cemento entre los años 2012 y 2019, lo que indica una expansión en la producción.Sin embargo, se evidencia una desaceleración en el año 2020 ocasionado por la pandemia. Adicionalmente, se observa que los despachos de cemento luego de la desaceleración del 2020 tuvieron un periodo de expansión. A partir de los años 2023 - 2024 se presenta un declive, lo que sugiere una desaceleración en la producción de cemento.

Extraer la tendencia DECEM

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
DECEM_vec <- as.numeric(DECEM_ts)
tendencia_DECEM <- as.numeric(stl_decomp_DECEM$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_DECEM <- ggplot() +
  geom_line(aes(x = fechas, y = DECEM_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_DECEM, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Gráfico 8 DECEM: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Toneladas") +
  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_DECEM)

Como se puede analizar en el gráfico 8, al comparar la serie original con la tendencia de los despachos de cemento se puede observar una expansión constante hasta 2019, con un aumento progresivo en la distribución y exportación.

La demanda de cemento aumentó hasta el año 2019, sin embargo en los últimos años la tendencia es descendente, lo que puede estar relacionada con una desaceleración del sector construcción o reducción del gasto en infraestructura.

Extraer la tendencia CONCRETO

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
CONCRETO_vec <- as.numeric(CONCRETO_ts)
tendencia_CONCRETO <- as.numeric(stl_decomp_CONCRETO$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_CONCRETO <- ggplot() +
  geom_line(aes(x = fechas, y = CONCRETO_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_CONCRETO, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
  ggtitle("Gráfico 9. CONCRETO: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Metros cúbicos") +
  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_CONCRETO)

En el gráfico 9 se puede observar en la tendencia un crecimiento inicial sólido hasta los años 2015-2016, seguido de un descenso progresivo hasta 2020 y una fuerte caída durante la pandemia.Posteriormente, se identifica una recuperación post-pandemia, aunque sin alcanzar los niveles máximos previos con una posible desaceleración reciente, lo que podría reflejar ajustes en la demanda de concreto.

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

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_PNCEM <- (PNCEM_ts[(13:length(PNCEM_ts))] / PNCEM_ts[1:(length(PNCEM_ts) - 12)] - 1) * 100
tasa_tendencia_PNCEM <- (tendencia_PNCEM[(13:length(tendencia_PNCEM))] / tendencia_PNCEM[1:(length(tendencia_PNCEM) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_PNCEM))
## [1] 144
print(length(tasa_crecimiento_PNCEM))
## [1] 144
print(length(tasa_tendencia_PNCEM))
## [1] 144
# Gráfico de la tasa de crecimiento anual
grafico_crecimiento_PNCEM <- ggplot() +
  geom_line(aes(x = fechas_corregidas_PNCEM, y = tasa_crecimiento_PNCEM), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_PNCEM, y = tasa_tendencia_PNCEM), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Gráfico 10 producción de cemento: Tasa de crecimiento anual % de la serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_PNCEM)

En el gráfico 10, la tasa de crecimiento de la tendencia se mantiene cercana a 0%, lo que indica un crecimiento estable o con pequeñas variaciones y algunas fluctuaciones menores, pero sin cambios drásticos. En el 2020 Se observa una caída a valores negativos, reflejando una contracción en la producción de CEMEX y del sector en general, lo que concuerda con la crisis global por la pandemia.

En el año 2021 se identifica una recuperación tras la caída del año 2020, hay un pico de crecimiento, con un aumento del 18.24%. Esto puede deberse a la reactivación de la construcción tras el impacto inicial del 2020.

No obstante, la tasa de crecimiento vuelve a estabilizarse cerca del 0%, pero en los últimos años muestra una leve desaceleración del -5.91%

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_DECEM <- (DECEM_ts[(13:length(DECEM_ts))] / DECEM_ts[1:(length(DECEM_ts) - 12)] - 1) * 100
tasa_tendencia_DECEM <- (tendencia_DECEM[(13:length(tendencia_DECEM))] / tendencia_DECEM[1:(length(tendencia_DECEM) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_DECEM))
## [1] 144
print(length(tasa_crecimiento_DECEM))
## [1] 144
print(length(tasa_tendencia_DECEM))
## [1] 144
# Gráfico de la tasa de crecimiento anual
grafico_crecimiento_DECEM <- ggplot() +
  geom_line(aes(x = fechas_corregidas_DECEM, y = tasa_crecimiento_DECEM), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_DECEM, y = tasa_tendencia_DECEM), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Gráfico 11 despachos de cemento: Tasa de crecimiento anual % de la serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_DECEM)

Como se puede observar en el gráfico 11 la tasa de crecimiento anual de los despachos de cemento está en una fase estable y predecible antes del año 2020.Sin embargo, es clave analizar una tendencia decreciente en los últimos años, visualizándose una contracción que llega a -6.21%.

Estos datos reflejan una disminución sostenida en los despachos de cemento en Colombia en los últimos años, lo que podría estar relacionado con una desaceleración en el sector de la construcción y una menor demanda de materiales.

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_CONCRETO <- (CONCRETO_ts[(13:length(CONCRETO_ts))] / CONCRETO_ts[1:(length(CONCRETO_ts) - 12)] - 1) * 100
tasa_tendencia_CONCRETO <- (tendencia_CONCRETO[(13:length(tendencia_CONCRETO))] / tendencia_CONCRETO[1:(length(tendencia_CONCRETO) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_CONCRETO))
## [1] 144
print(length(tasa_crecimiento_CONCRETO))
## [1] 144
print(length(tasa_tendencia_CONCRETO))
## [1] 144
# Gráfico de la tasa de crecimiento anual
grafico_crecimiento_CONCRETO <- ggplot() +
  geom_line(aes(x = fechas_corregidas_CONCRETO, y = tasa_crecimiento_CONCRETO), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_CONCRETO, y = tasa_tendencia_CONCRETO), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Gráfico 12 producción de concreto premezclado: Tasa de crecimiento anual % de la serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_CONCRETO)

En el gráfico 12 se refleja que el crecimiento anual de la tendencia en la producción de concreto premezclado en Colombia no es muy variable a lo largo del tiempo.También se puede evidenciar una caída de la tasa de crecimiento en 2020 por pandemia COVID-19, pero es seguida de una recuperación rápida entre los años 2021 - 2022, lo que indica la resiliencia del sector.

La producción de concreto premezclado ha experimentado una tendencia decreciente en los últimos años 2023 - 2024, alcanzando los -5.20%, influenciado por una menor demanda en sectores clave como la vivienda y las obras civiles.

Modelo Arima

El modelo ARIMA (Autoregressive Integrated Moving Average) es una de las metodologías más utilizadas en la predicción de series temporales. Su capacidad para modelar datos secuenciales permite analizar tendencias, patrones estacionales y fluctuaciones en el tiempo.

En el caso PNCEM (Producción de Cemento) en Colombia, esta variable está influenciada por factores como la demanda del sector de la construcción, la inversión en infraestructura, las políticas gubernamentales y las condiciones macroeconómicas. Al aplicar un modelo ARIMA, podemos identificar patrones históricos y generar pronósticos que ayuden en la planificación de la producción y la toma de decisiones estratégicas.

Metodología Box-Jenkins para aplicar un modelo ARIMA

La metodología Box-Jenkins es un enfoque sistemático para construir modelos ARIMA con el objetivo de analizar y pronosticar series de tiempo. Fue desarrollada por George Box y Gwilym Jenkins y se basa en cuatro etapas clave:

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

Antes de empezar a aplicar la metodlogía BOX-JENKINS a la variable PNCEM,se debe dividir el conjunto de datos de prueba y entrenamiento

💡 Esta división es clave en modelos predictivos para evitar sobreajuste y evaluar el rendimiento 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: octubre 2024-diciembre 2024 

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

Paso 1: Identificación del modelo

Identificar estacionariedad

En el análisis de series de tiempo, una serie es estacionaria si su comportamiento es constante a lo largo del tiempo, es decir:

✅ Su promedio no cambia con el tiempo. ✅ Su variabilidad (qué tanto fluctúa) se mantiene estable. ✅ Su relación con valores pasados es siempre la misma.

¿Cómo saber si una serie es estacionaria?

1️⃣ Observando un gráfico Si el gráfico de la serie muestra una tendencia creciente o decreciente, o si las variaciones se hacen más grandes con el tiempo, la serie probablemente no es estacionaria.

2️⃣ Usando la prueba de Dickey-Fuller Aumentada (ADF) Es una prueba estadística que nos dice si la serie tiene una tendencia fuerte. Si el p-valor de la prueba es mayor a 0.05, significa que la serie no es estacionaria.

Test de Dickey-Fuller

El test de Dickey-Fuller aumentado (ADF) se usa para verificar si una serie temporal es estacionaria, es decir, si sus propiedades estadísticas (media y varianza) permanecen constantes en el tiempo.

HO: Serie no estacionaria HI: Serie estacionaria

¿Qué significa el p-valor?

Si el p-valor es bajo (< 0.05) → Rechazamos la hipótesis nula y concluimos que la serie es estacionaria. Si el p-valor es alto (> 0.05) → No podemos rechazar la hipótesis nula, lo que indica que la serie no es estacionaria.

A continuación se aplica el test ADF para validar estacionariedad en el conjunto de entrenamiento de la variable escogida PNCEM (producción de cemento), que es la elegida para pronosticar:

library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable PNCEM (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_ts
## Dickey-Fuller = -4.5156, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

A pesar de que el resultado del augmented Dickey-Fuller (ADF) arrojó un valor P-value por debajo de 0.05 indicando que la serie es estacionaria, se realizará diferenciación de niveles debido a que en el gráfico se evidencia que no es estacionaria, por lo cual se hace necesario aplicar una diferenciación.

Diferenciación en niveles variable PNCEM

#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable PNCEM , una sola vez:
train_diff <- diff(train_ts, differences = 1) 

A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable PNCEM para ver graficamente el cambio o ajuste:

# Graficar la serie original en un gráfico separado
  p2 <- ggplot(data.frame(Tiempo = time(train_ts), PNCEM = as.numeric(train_ts)), aes(x = Tiempo, y = PNCEM)) +
    geom_line(color = "blue") +
    ggtitle("Gráfico 13 PNCEM:Serie Original") +
    xlab("Tiempo") + ylab("Toneladas de cemento")
  
  ggplotly(p2)  # Convertir en gráfico interactivo

En este gráfico se puede observar las toneladas de cemento producidas a lo largo del tiempo en un periodo comprendido del 01 de enero del 2012 al 01 de diciembre del 2024 en su versión original.

# Graficar la serie diferenciada en un gráfico separado
  p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], PNCEM_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = PNCEM_Diff)) +
    geom_line(color = "red") +
    ggtitle("Gráfico 14 PNCEM:Serie Estacionaria (Una diferenciación)") +
    xlab("Tiempo") + ylab("variable PNCEM diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

En este gráfico de la variable PNCEM con la difrenciación aplicada se logra evidenciar que se suavizan los datos en comparación con el gráfico 13 de la versión original, pero aún no se logra que la serie sea estacionaria, aún se observan muchas fluctuaciones.

Diferenciación en logaritmo

La aplicación de logaritmos en series de tiempo se realiza principalmente para lograr que la serie sea más estable, lineal y estacionaria. Esta transformación es relevante porque permite modelar mejor las series que siguen un crecimiento exponencial y facilita la aplicación de técnicas estadísticas que requieren estacionariedad.

A continuación se aplica la diferenciación logarítimica y la varible PNCEM ahora se llama train_diff_log:

# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación

train_diff_log <- diff(log(train_ts), differences = 1)

Ahora graficamos la serie orignal versus la serie diferenciada una vez con logaritmo

# Graficar la serie original en un gráfico separado
  p2 <- ggplot(data.frame(Tiempo = time(train_ts), PNCEM = as.numeric(train_ts)), aes(x = Tiempo, y = PNCEM)) +
    geom_line(color = "blue") +
    ggtitle("Gráfico 15 PNCEM:Serie Original") +
    xlab("Tiempo") + ylab("PNCEM")
  
  ggplotly(p2)  # Convertir en gráfico interactivo

En este gráfico nuevamente se puede observar las toneladas de cemento producidas a lo largo del tiempo en un periodo comprendido del 01 de enero del 2012 al 01 de enero del 2024 en su versión original.

# Graficar la serie diferenciada en un gráfico separado
  
  p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], micro_Diff = as.numeric(train_diff_log)), aes(x = Tiempo, y = micro_Diff)) +
    geom_line(color = "red") +
    ggtitle("Gráfico 16 PNCEM:Serie Estacionaria (Una diferenciación en logaritmo)") +
    xlab("Tiempo") + ylab("PNCEM diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

En este gráfico se aplico la diferenciación logarítmica para lograr una mejor serie estacionaria para la variable PNCEM, pero se puede observar que aún existen fluctuaciones en los datos, aunque en comparación con el gráfico 14 en el cual se le aplico una diferenciación de niveles, se puede concluir que la diferenciación logarítmica muestra mejores resultados.

Ahora se prueba estacionariedad en la serie diferenciada (en nivel y logaritmo)

# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
  adf_test_diff <- adf.test(train_diff)
  print(adf_test_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff
## Dickey-Fuller = -7.4517, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

A pesar de que el resultado del augmented Dickey-Fuller (ADF) arrojó un valor p-value por debajo de 0.05 indicando que la serie es estacionaria, el gráfico 14 muestra que la serie no es estacionaria.

# Segunda prueba de estacionariedad sobre la serie diferenciada en logaritmo
  adf_test_diff_log <- adf.test(train_diff_log)
  print(adf_test_diff_log)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff_log
## Dickey-Fuller = -7.6005, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Lo mismo ocurre con este resultado del augmented Dickey-Fuller (ADF) arrojó un valor p-value por debajo de 0.05 indicando que la serie es estacionaria, como se observa en el gráfico 16, se logró suavizar la serie en los primeros y últimos años del periodo. Sin embargo, la diferenciación logarítmica no logró estabilizar los años 2020- 2021.

El p-value ya es menor a 0.05 con una primera diferencia de niveles y logaritmica, Por tanto el valor que puede tomar d es igual a 1.

Identificación manual de p y q

ACF (Autocorrelation Function)

Muestra la correlación de la serie con sus rezagos.

Ayuda a determinar el parámetro q en un modelo ARIMA(p, d, q).

PACF (Partial Autocorrelation Function)

Muestra la correlación parcial entre la serie y un rezago específico, eliminando el efecto de rezagos intermedios.

Ayuda a determinar el parámetro p en un modelo ARIMA(p, d, q).

library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff_log, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")

ggplotly(acf_plot)
ggplotly(pacf_plot)

Interpretación correlogramas

Se puede observar que los valores que podrian tomar p y q serian:

p=1*

q=1* q2 q=3 q=4

(P óptimo=1) (q óptimo=1)

El modelo óptimo para esta variable seria (1,1,1)

Paso 2. Estimación manual del modelo

Estimación del modelo identificado (1,1,1)

# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(1,1,1))
summary(manual_arima_model)
## 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

El summary del modelo Arima manual de la variable producción de cemento indica que la aplicación de este presenta un sesgo promedio subestimado de 10.449 toneladas, reflejando un sesgo sistemático considerable de manera consistente. Adicionalmente, se analizaron otras métricas para obtener una evaluación completa del rendimiento del modelo, tales como: RMSE, MAE, MPE, MAPE, ACF1, estas métricas se encuentran alineadas en los resultados:

RMSE indica 105.259 errores en el pronóstico del modelo.

MAE: el error promedio absoluto entre las predicciones y los valores observados es de 72.659

MAPE: 8,9% de errores en el modelo.

MPE: muestra un valor negativo de -1,5% lo que refleja que el modelo sobreestima los valores reales (las predicciones son mayores que los valores observados).

ACF1: indica que los valores consecutivos de la serie temporal están casi relacionados entre sí, con una correlación débilmente negativa.

De acuerdo a lo anterior, el modelo Arima manual no es aceptable para esta variable, por lo que requiere ajustes adicionales tales como: un modelo diferente o diferenciación logarítmica.

Significancia de coefcientes

library(lmtest)

# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)
## 
## 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

Interpretación significancia coeficientes

ar1, es altamente significativo (***), lo que significa que este coeficiente tiene un impacto importante en el modelo.

 checkresiduals(manual_arima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 45.196, df = 22, p-value = 0.002506
## 
## Model df: 2.   Total lags used: 24

Paso 3. Validación de residuos del modelo estimado manual

La validación de residuos es crucial para determinar si el modelo ARIMA es adecuado o si necesita mejoras. El objetivo es verificar que los residuos (errores de predicción) se comporten como ruido blanco, es decir, sin patrones detectables.

1. Serie de residuos (gráfico superior)

Muestra cómo se comportan los errores a lo largo del tiempo.

Idealmente, deberían oscilar alrededor de cero sin tendencias evidentes ni grandes acumulaciones de error.

Problema posible: Se observan picos alrededor del 2020, lo que podría indicar cambios estructurales en los datos (como efectos de la pandemia en la producción de cemento).

Función de Autocorrelación (gráfico inferior izquierdo, ACF de residuos)

Si el modelo es adecuado, los residuos no deben mostrar correlaciones significativas en el tiempo.

Interpretación:

La mayoría de las barras están dentro de las líneas azules (intervalos de confianza). Sin embargo, algunos rezagos parecen salir del rango, lo que sugiere que aún puede haber estructura no capturada en los datos.Esto indica que el modelo podría mejorarse.

Histograma de residuos con ajuste normal (gráfico inferior derecho)

Sirve para verificar si los errores siguen una distribución normal, lo cual es un supuesto clave en ARIMA.

Interpretación:

La curva roja representa la distribución normal teórica. Los residuos se acercan a la normalidad, pero hay algunos valores extremos (colas más gruesas de lo esperado).Esto indica que puede haber eventos atípicos o datos no bien explicados por el modelo en este caso pandemia.

Conclusión y acciones recomendadas

✅ El modelo ARIMA(1,1,1)parece razonablemente bueno, pero tiene algunas señales de que podría mejorarse.

⚠️ Posibles mejoras:

Revisar la estructura del modelo: Se pueden probar otros órdenes ARIMA o incluso modelos más complejos como SARIMA o modelos con variables exógenas (ARIMAX).

Manejo de valores atípicos: Considerar incluir un término de intervención si eventos como la pandemia afectaron la serie.

Transformación de datos: Si los residuos no son normales, una transformación logarítmica o Box-Cox puede ayudar.

Recordar: El supuesto de normalidad significa que los errores o residuos de un modelo deben seguir una distribución normal (o “campana de Gauss”). Si los errores son normales, podemos hacer predicciones más confiables y usar ciertas pruebas estadísticas que asumen esta propiedad.

Paso 4. Pronóstico (modelo manual)

El modelo sigue la tendencia general, pero tiene un sesgo de sobreestimación.

Si la serie tiene patrones estacionales fuertes y no están bien capturados, el modelo puede fallar en prever las fluctuaciones.

Pronóstico en el test de prueba (oct, nov y dic 2024) y gráfico

#Aquí se crea el pronóstico con el modelo ARIMA manual y se guarda en una nueva riable u objeto "manual_forecast"
manual_forecast <- forecast(manual_arima_model, h = length(test_ts)) #Se generan tantos pronósticos como valores tenga el conjunto de prueba (test_ts).

# Crear dataframe para gráfico interactivo del pronóstico manual
manual_forecast_data <- data.frame(Tiempo = time(manual_forecast$mean), ## Extrae las fechas del pronóstico
                                   Pronostico = as.numeric(manual_forecast$mean), ## Valores pronosticados
                                   Observado = as.numeric(test_ts)) ## Valores reales de la serie

# Graficar pronóstico manual junto con los valores observados reales
p_manual <- ggplot(manual_forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico Manual")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Gráfico 17 PNCEM:Pronóstico Manual vs Observado") +
  xlab("Tiempo") + ylab("PNCEM")

ggplotly(p_manual)

De acuerdo al gráfico 17, el pronóstico manual no anticipó la caída fuerte que se observó en los datos reales.

El valor final del PNCEM observado y el pronosticado son similares, lo que sugiere que, a pesar de la variabilidad en los valores reales, la tendencia general del pronóstico fue acertada.

El pronóstico manual fue razonable en tendencia, pero no capturó la variabilidad en los datos observados.

Métricas de evaluación del modelo manual dentro del periodo de prueba (oct,nov y dic 2024)

# Calcular métricas de evaluación del modelo manual
mae_manual <- mean(abs(manual_forecast$mean - test_ts), na.rm = TRUE)
rmse_manual <- sqrt(mean((manual_forecast$mean - test_ts)^2, na.rm = TRUE))

# Mostrar métricas de evaluación del modelo manual
cat("MAE Manual: ", mae_manual, "\n")
## MAE Manual:  23001.89
cat("RMSE Manual: ", rmse_manual, "\n")
## RMSE Manual:  29742.55

MAE (Mean Absolute Error)

Indica el error promedio en unidades de la variable, es decir, que en promedio el modelo se equivoca en 23.002 toneladas de cemento al realizar su predicción.

RMSE (Root Mean Squared Error)

Penaliza más los errores grandes debido a la elevación al cuadrado antes de calcular la raíz. Su valor de 29.743 toneladas de cemento, lo que sugiere que hay algunos errores grandes que afectan la precisión del modelo.

A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic 2024

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

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

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

# Mostrar la tabla
print(forecast_table_manual)
##     Tiempo Observado Pronosticado
## 1 2024.750   1180815      1133340
## 2 2024.833   1121219      1141154
## 3 2024.917   1145357      1143761

Como se puede observar en este resultado la presición del modelo necesita mejoras, dado que en el primer periodo al calcular la diferencia entre lo observado y lo pronosticado se evidencia que el modelo subestimó la demanda en -47.475 toneladas de producción de cemento, por lo cual se deben implementar mejoras en este modelo para ser más acertivos en el pronóstico.

Ahora pronosticamos fuera del periodo de análisis

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_manual <- forecast(manual_arima_model, h = length(test_ts) + 1)

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

# Mostrar el pronóstico completo
print(next_month_forecast_manual)
##     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_manual, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2025 = 1144630.62871857"

Otra forma para calcular un valor futuro (fuera de muestra)-Modelo manual

Si no hay test y solo quieres el siguiente punto u observación, el código a usar seria algo asi:

# Pronosticar octubre 2024
future_forecast_manual <- forecast(manual_arima_model, h = 1)

# Extraer el valor específico de octubre 2024
forecast_oct2024 <- future_forecast_manual$mean[1]
print(paste("Pronóstico para oct 2024:", forecast_oct2024))
## [1] "Pronóstico para oct 2024: 1133339.64973436"

Modelo Arima automático

Usa la función auto.arima() de forecast en R para seleccionar automáticamente los mejores parámetros (p,d,q).

✅ Ventajas:

✔ Optimización automática: Encuentra los valores óptimos de ARIMA sin intervención manual. ✔ Ahorra tiempo: Útil cuando hay muchas series a modelar. ✔ Evita sesgo humano: Reduce el riesgo de elegir un modelo incorrecto por falta de experiencia. ✔ Incluye corrección por estacionalidad si se usa con seasonal = TRUE. ✔ Suele funcionar bien en la mayoría de los casos, ya que usa criterios como AIC/BIC para optimizar.

❌ Desventajas: ❌ Puede no ser el mejor modelo posible, ya que depende del criterio de selección. ❌ Menos interpretabilidad: No siempre es claro por qué eligió ciertos parámetros. ❌ Puede ignorar conocimiento experto sobre la serie o factores externos.

Identificación automática del modelo ARIMA

library(forecast)

# Ajustar un modelo ARIMA automático sin estacionalidad
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

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

ar1, es altamente significativo (***), lo que significa que este coeficiente tiene un impacto importante en el modelo.

# Ajuste del modelo ARIMA(1,1,1) sin parte estacional
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.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

El summary del modelo Arima automático de la variable producción de cemento nos indica que la aplicación de este presenta un sesgo sistemático dado que el resultado de la diferencia promedio entre los valores reales y predichos están en 10.449 de error.Adicionalmente, se analizaron otras métricas para obtener una evaluación completa del rendimiento del modelo automático, tales como: RMSE, MAE, MPE, MAPE, ACF1.

RMSE indica 105.259 errores en el pronóstico del modelo.

MAE: el error promedio absoluto entre las predicciones y los valores observados es de 72.659

MAPE: 8,9% de errores en el modelo.

MPE: muestra un valor negativo de -1,5% lo que refleja que el modelo sobreestima los valores reales (las predicciones son mayores que los valores observados).

ACF1: indica que los valores consecutivos de la serie temporal están casi relacionados entre sí, con una correlación débilmente negativa.

De acuerdo a lo anterior, el modelo Arima automático no es aceptable para esta variable, por lo que requiere ajustes adicionales tales como: un modelo diferente o diferenciación logarítmica.

En conclusión, el modelo arima automático arrojó los mismos resultados al modelo arima manual.Por ende este modelo tampoco es óptimo para realizar el pronóstico de la variable de producción de cemento.

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

Pronóstico modelo ARIMA automático

# 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("Gráfico 18 Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("Microempresas")

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

Como se puede observar en el gráfico 18, el pronóstico automático vs lo observado muestra el mismo resultado de la aplicación del modelo manual. Es decir, que el modelo no capturó la variabilidad de los datos, lo que indica que el modelo necesita ajustes dado que no está prediciendo de forma precisa.

# 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

Los resultados de este pronóstico es igual a los resultados del modelo manual.

Ahora pronosticamos fuera del periodo de análisis

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

Otra forma para calcular un valor futuro (fuera de muestra)

# Pronosticar el mes de octubre 2024
future_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = 1)

# Extraer el valor específico de octubre 2024
forecast_oct2024_auto <- future_forecast_auto$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_auto))
## [1] "Pronóstico para octubre 2024: 1133339.64973436"

Identificación automática del modelo Sarima

Un SARIMA (Seasonal ARIMA) es un modelo que combina:

✅ Un modelo ARIMA tradicional para capturar relaciones entre valores pasados y errores. ✅ Un componente estacional, útil cuando los datos muestran patrones repetitivos en el tiempo (por ejemplo, ventas trimestrales, datos mensuales de temperatura, etc.).

Básicamente, es como un ARIMA mejorado que puede manejar ciclos estacionales.

# Identificación automática del mejor modelo ARIMA
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

De acuerdo a los parámetros arrojados en este estudio del modelo Sarima, se puede observar que:

El coeficiente AR(1) indica un valor positivo (0.3748) lo que muestra una correlación positiva entre los valores de la serie temporal en los períodos consecutivos.

El error estándar de 0,0966 es relativamente pequeño en comparación con el coeficiente, lo que indica que el valor estimado para ar1 es bastante confiable.

Adicionalmente, el error estándar del coeficiente MA es de 0,0477, que también es un valor relativamente pequeño y presenta una estimación confiable del parámetro.

El coeficiente de SMA refleja un valor positivo de 0,3120 indica que el error estacional en el período anterior (12 períodos antes) tiene una influencia positiva sobre el valor actual; y su error estándar de 0,00782 también es razonablemente bajo.

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(1,1,1)(0,0,1)[12]
darima <- Arima(train_ts, 
                order = c(1, 1, 1),  # (p,d,q) -> (1,1,1)
                seasonal = list(order = c(0, 0, 1),  # (P,D,Q) -> (0,0,1)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(1,1,1)(0,0,1)[12] 
## 
## Coefficients:
##          ar1      ma1    sma1
##       0.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

El modelo Sarima se descompone en dos partes:

Parte no estacional

AR(1) = 0.3748 Indica que el valor actual de la serie está influenciado en un 37.48% por su valor anterior.

MA(1) = -0.9184 Indica que los errores pasados tienen un peso importante en la predicción, con un coeficiente negativo alto, lo que sugiere que los errores grandes en un período tienden a compensarse en el siguiente.

parte estacional 12 meses

SMA(1) = 0.3120 Indica que hay un patrón de media móvil en la estacionalidad anual, con un peso del 31.2% en los errores estacionales pasados.

A continuación mencionamos las métricas de error del modelo más significativas.

ME (Error Medio) = 9371 Diferencia promedio entre el valor real y el pronosticado. como este valor es lejano a 0 indica que hay sesgo.

MPE (Error porcentual medio) = -1.49% Indica que hay una ligera subestimación en promedio.

MAPE (Error absoluto porcentual medio) = Un error del 8.36% en promedio. Un MAPE inferior al 10% es un modelo relativamente bueno para pronóstico.

En conclusión, el modelo captura la tendencia y la estacionalidad correctamente, aunque tiene cierto margen de error, pero es óptimo para ponóstico.

También se evidencia que los coeficientes indican que los valores pasados y los errores tienen un peso significativo en la predicción.

# 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

Pronóstico con el modelo SARIMA

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

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

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

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

El pronóstico comienza con un valor más bajo que el observado, pero ambos muestran una tendencia a la baja.

Aunque ambos alcanzan su punto más bajo en el mismo período de tiempo, el valor mínimo observado es ligeramente superior al pronosticado.

El pronóstico predice una recuperación más rápida y agresiva que la realidad observada.

En conclusión el modelo Sarima es aceptable para realizar un pronóstico de la variable PNCEM (producción de cemento).

# 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

El modelo SARIMA está funcionando razonablemente bien, aunque en algunos casos subestima o sobreestima los valores reales.

Ahora pronosticamos fuera del periodo de análisis

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

# 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 2024.750    1147233
## 2 2024.833    1117421
## 3 2024.917    1171977
## 4 2025.000    1101935
# Extraer solo el valor del mes 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"

El pronóstico muestra pequeñas fluctuaciones, con una caída prevista a inicios de 2025.

Otra forma para calcular un valor futuro (fuera de muestra)

# Identificación automática del mejor modelo ARIMA
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

El ma1 es muy significativo (error pequeño), lo que indica que la parte de media móvil tiene un gran impacto.

El sma1 también es relevante (error moderado).

ar1 tiene más incertidumbre, pero sigue siendo significativo.

# Pronosticar octubre 2024
future_forecast_sarima <- forecast(auto_arima_model, h = 1)

# Extraer el valor específico de octubre 2024
forecast_oct2024_sarima <- future_forecast_sarima$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_sarima))
## [1] "Pronóstico para octubre 2024: 1147233.34511743"

El modelo SARIMA ha predicho que el valor para octubre de 2024 será 1.147.233.

Hallazgos clave para la compañía

Crecimiento y Desaceleración del Sector construcción en el periodo (2012-2024):

En los años 2012-2016 hubo un Crecimiento sostenido impulsado por programas de vivienda social como Mi Casa Ya y proyectos de infraestructura 4G.

La construcción en Colombia entre 2017 y 2019 enfrentó un periodo de desaceleración debido a la reducción de la inversión en vivienda e infraestructura, problemas de financiamiento y aumento de costos. Sin embargo, hacia el final del período, algunas medidas gubernamentales comenzaron a generar un repunte en la actividad.

En el 2020-2021 se presentó una fuerte caída debido a la pandemia y el paro nacional (cierres de obras, problemas con el transporte de materiales, aumento de costos de materiales, entre otros).

Del año 2022 al 2024 se obtuvo una recuperación parcial en obras civiles, pero dificultades en vivienda y edificaciones privadas tales como el aumento de la inflación que elevó los costos y afectó la demanda de cemento.

En conclusión, entre 2012 y 2024 la producción de cemento en Colombia mostró un crecimiento inicial significativo, seguido de una expansión impulsado por programas de vivienda. Sin embargo, en los años recientes, factores económicos y una disminución en la actividad constructiva condujeron a una desaceleración en la producción y los despachos de cemento, evidenciando la necesidad de estrategias para revitalizar el sector.

Impactos del pronóstico en el sector construcción

De acuerdo a los resultados obtenidos del pronóstico, Se observa un comportamiento fluctuante con descensos en noviembre 2024 y enero 2025, lo que muestra una desaceleración en la demanda de materiales de construcción que podría indicar lo siguiente:

-Reducción en la actividad constructora, especialmente en proyectos de vivienda y edificaciones comerciales.

-Retrasos o cancelaciones en la construcción de nuevas obras debido a la menor disponibilidad de materiales o aumento de costos.

-Menor inversión en infraestructura, si la demanda de cemento se mantiene baja.

-Aumento en costos de materiales, si la oferta disminuye y la demanda se mantiene.

-Reducción de la contratación de trabajadores en el sector.

En conclusión la construcción es un motor clave en la economía colombiana. Una desaceleración del sector podría afectar el PIB y ralentizar el crecimiento de otras industrias relacionadas (acero, vidrio, cerámica).

Recomendaciones estratégicas CEMEX Colombia S.A.

Ante el pronóstico de una desaceleración en la producción de cemento y posible impacto en el sector construcción en 2025, la compañía CEMEX COLOMBIA S.A.debe adaptarse estratégicamente para mitigar riesgos y mantener la rentabilidad.

Estrategias Financieras y de Costos

-Optimización del Presupuesto tales como reducir gastos innecesarios en la operación y mejorar la eficiencia de costos.

-Negociar precios con proveedores y asegurar contratos de suministro a largo plazo.

-Diversificar fuentes de financiamiento para evitar dependencia de créditos bancarios con tasas altas.

Diversificación de Ingresos

-Explorar nuevas líneas de negocio como remodelación y mantenimiento.

-Expandirse hacia sectores menos afectados como infraestructura pública o proyectos de energía renovable.

Gestión de Riesgo Financiero

-Monitorear indicadores económicos y ajustar planes de inversión.

-Mantener reservas de liquidez para hacer frente a posibles caídas en la demanda.

Innovación y Eficiencia Operativa

  • Uso de Nuevas Tecnologías y Construcción Sostenible.

Automatización y Digitalización

-Aplicar BIM (Building Information Modeling) para mejorar la planificación y optimizar recursos.

-Usar herramientas de gestión de proyectos en la nube para mayor eficiencia operativa.

Estrategias Comerciales y de Mercado

-Segmentación de Clientes y Nuevas Oportunidades.

-Enfocarse en nichos de mercado en crecimiento como vivienda de interés social (VIS). Explorar alianzas con gobiernos locales para participar en proyectos de infraestructura pública.

-Expandir operaciones hacia ciudades secundarias donde la demanda aún se mantenga estable.

Estrategias de Fidelización y Marketing

-Ofrecer planes de financiamiento flexibles a clientes para facilitar ventas.

-Implementar estrategias de marketing digital y ventas en línea para llegar a más compradores.

En conclusión, Si bien la desaceleración en 2025 puede afectar la industria, las empresas que se adapten con innovación, diversificación y control de costos tendrán mayores posibilidades de mantenerse competitivas, por tal motivo CEMEX deberá implementar las estrategias anteriormente descritas con el fin de fortalecer su posición en el mercado, optimizar sus operaciones y prepararse para una eventual reactivación del sector. Al adoptar un enfoque proactivo e innovador, CEMEX podrá mitigar los efectos de la desaceleración, garantizar su sostenibilidad a largo plazo y aprovechar nuevas oportunidades de crecimiento cuando el mercado se recupere.