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.

Introducción El análisis de la producción de pollo, la producción de huevo y la importación de cereales en Colombia entre 2012 y 2024 permite comprender la dinámica integral de la industria avícola, uno de los sectores más relevantes para la seguridad alimentaria del país. Estas tres variables están estrechamente relacionadas: la producción de pollo y huevo depende en gran medida del suministro y costo de cereales como maíz y soya, los cuales son importados en un alto porcentaje para la formulación de alimentos balanceados.

Cargar base de datos

library(readxl)
data_col <- read_excel("~/PERSONAL/MAESTRIA/CASO 2/Base Caso2-.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric"))

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

Variable 1

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

Variable 2

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

Variable 3

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

Extracción de señales

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

library(ggplot2)
library(plotly)

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

# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)), 
                                      y = variable1)) +
  geom_line(color = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Serie original: POLLO") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Extracción señales variable 1 Tendencia alcista en los últimos años.

# 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: POLLO",
       x = "Tiempo",
       y = "Valor")

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

Análisis a la descomposición temporal de la producción de pollo en Colombia 2012-2024 - en toneladas Estacionalidad (Línea roja) Patrón repetitivo anual: Cada año se aprecia un pico y un valle claros, reflejando estacionalidad pronunciada. Picos: Posiblemente en junio-julio y noviembre-diciembre, coincidiendo con temporadas de mayor consumo (vacaciones y fin de año). Valles: En enero-febrero y a veces en mitad de año, típico del ciclo de demanda y producción.

Componente residual (Línea verde) Mayor dispersión en 2020: Hay valores atípicos negativos muy marcados en 2020, que reflejan choques imprevistos (pandemia, restricciones logísticas, costos de insumos). Después de 2021, los residuos vuelven a ser relativamente estables y centrados en torno a cero, señal de recuperación y normalización.

Serie original (Línea azul) Combina los efectos anteriores: crecimiento a largo plazo, estacionalidad anual y perturbaciones puntuales. Claramente hay un valle profundo en 2020, que concuerda con el residuo negativo y la pausa en la tendencia.

Tendencia (Línea morada) Crecimiento sostenido: Desde 2012 hasta 2019 la tendencia es ascendente constante, pasando aproximadamente de 92.688 a cerca de 141.000 toneladas mensuales. Impacto en 2020: Se observa un quiebre leve o desaceleración alrededor de 2020 (coincidente con la pandemia y disrupciones en la cadena de suministros). Recuperación posterior: Entre 2021 y 2023 retoma la tendencia positiva, pero con un ritmo más moderado, llegando a estabilizarse ligeramente en 2024.

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: iMP_CEREALES",
       x = "Tiempo",
       y = "Valor")

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

Análisis a la descomposición temporal de las importaciones de cereales en Colombia 2012-2024 -en miles de dólares Estacionalidad (Línea roja) Patrón repetitivo marcado: Ciclo anual con picos y valles claros, probablemente relacionados con programación de compras e importaciones anticipadas antes de temporadas de alta demanda (alimentación animal). Picos frecuentes posiblemente en el primer semestre (cuando se acumulan inventarios para mitad de año).

Residuos (Línea verde) Alta variabilidad durante todo el período, con valores positivos y negativos amplios. 2021-2022 muestra desviaciones más pronunciadas, lo que coincide con choques globales (crisis logística, incremento en precios internacionales, pandemia). Serie Original (Línea azul) Refleja claramente la volatilidad de importaciones: Subidas y bajadas frecuentes, pero con una clara tendencia alcista hasta 2022. Posterior caída hacia 2024, lo que sugiere una corrección.

Tendencia (Línea morada) 2012-2018: Tendencia relativamente estable con ligeras fluctuaciones, oscilando entre 140.000 y 160.000. 2019-2022: Fuerte aumento, llegando a un pico alrededor de 2022 (~260.000). Esto coincide con: Aumento en los costos internacionales de granos. Dependencia del maíz y soya importados para la industria avícola. 2023-2024: Descenso moderado, estabilizándose cerca de 200.000, lo que indica ajuste tras el pico.

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: HUEVO",
       x = "Tiempo",
       y = "Valor")

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

Análisis a la descomposición temporal de la producción de huevo en Colombia 2012-2024 – en toneladas Estacionalidad (Línea roja) Presenta un patrón muy marcado y constante: Fluctuaciones regulares anuales con picos y valles predecibles. Esto indica que la producción responde a ciclos estacionales (posiblemente consumo en festividades y ajustes en la postura).

Residuos (Línea verde) Mayor variabilidad entre 2018-2022, con caídas pronunciadas en algunos periodos (posibles choques de precios de insumos y pandemia). Después de 2022, la variabilidad tiende a estabilizarse.

Serie Original (Línea azul) Refleja la combinación de tendencia ascendente con picos estacionales. El incremento post-2021 muestra que el sector logró resiliencia tras las disrupciones.

Tendencia (Línea morada) 2012-2016: Crecimiento sostenido, pasando de ~55.000 a 65.000. 2017-2019: Aceleración moderada, llegando a 70.000-75.000. 2020: Disminución ligera, coincidiendo con la pandemia y posibles problemas logísticos y de costos (alimento balanceado). 2021-2024: Recuperación con tendencia al alza, llegando a valores cercanos a 90.000-95.000. Conclusión: A pesar de choques puntuales, la tendencia es claramente creciente en el largo plazo, mostrando expansión sostenida del sector.

Analisis interaccion de las tres variables

Durante este periodo, se observa una tendencia de crecimiento sostenido tanto en la producción de pollo como en la de huevo, impulsada por el aumento del consumo interno y la competitividad del sector. Sin embargo, esta expansión ha estado condicionada por la volatilidad en las importaciones de cereales, especialmente ante choques externos como la pandemia de 2020, el aumento en los precios internacionales y las fluctuaciones en la tasa de cambio. El comportamiento estacional, presente en las tres variables, refleja la influencia de factores cíclicos en la oferta y la demanda, mientras que los choques residuales evidencian la vulnerabilidad del sistema productivo ante eventos globales y locales. Comprender estas interdependencias es fundamental para diseñar estrategias que garanticen la sostenibilidad, la estabilidad de precios y la seguridad alimentaria en el país.

# 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 La corrección es medianamente influyente

# 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("POLLO: Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Ton") +
  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)

Serie ajustada por estacionalidad medianamente suavizada y ajustada a original

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

Serie ajustada por estacionalidad medianamente suavizada y ajustada a original

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("HUEVO: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: POLLO 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)

Serie temporal de la producción de pollo (Variable 1) frente a su tendencia obtenida por una descomposición estacional. En los último años la tendencia es estable, con una ligera caida finalizando el año 2024.

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 debe comenzar desde la observación 13 ya que las anteriores no se pueden calcular y siempre debe cualcular 12 periodos atrás.

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

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

# Verificar longitudes
print(length(fechas_corregidas_var1))

## [1] 144

print(length(tasa_crecimiento_var1))

## [1] 144

print(length(tasa_tendencia_var1))

## [1] 144

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

library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable1:POLLO 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)

###Se rompe ciclo de crecimiento con el mínimo después de pandemia, llegando a -7% en el mes de Enero del año 2020. Después de esto se da un incremento finalizando el año 2021, se estabiliza cercano a números del 7% y cae finalizando el año 2022 debido a los bloqueos de vías que afectaron la logística y una fuerte caída en la inversión en activos biológicos (aves ponedoras) que redujo la oferta, sumado a la dinámica inflacionaria que afectó los costos de producción. En proyección sugiere una posible estabilización o caida de la tasa de crecimiento.

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

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

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

# Verificar longitudes
print(length(fechas_corregidas_var2))

## [1] 144

print(length(tasa_crecimiento_var2))

## [1] 144

print(length(tasa_tendencia_var2))

## [1] 144

# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable2 IMP_ CEREAL: 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)

La caida reflejada entre el año 2022 y 2023 debido a la volatilidad de los precios de los productos básicos, el aumento de la inflación, el fortalecimiento del dólar y el contexto de la guerra en Ucrania que afectó los precios y la oferta de alimentos a nivel mundial.

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

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_var3))

## [1] 144

print(length(tasa_crecimiento_var3))

## [1] 144

print(length(tasa_tendencia_var3))

## [1] 144

# Gráfico de la tasa de crecimiento anual variable 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 HUEVO: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

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

Modelo ARIMA

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

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

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

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

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

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

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

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

Identificación automática del modelo ARIMA

library(forecast)

# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)

# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)

## Series: train_ts 
## ARIMA(2,1,2) with drift 
## 
## Coefficients:
##           ar1     ar2      ma1      ma2     drift
##       -0.0113  0.3641  -0.1536  -0.7563  405.6258
## s.e.   0.1262  0.1159   0.0900   0.0871   71.5809
## 
## sigma^2 = 31352185:  log likelihood = -1525.54
## AIC=3063.09   AICc=3063.67   BIC=3081.23
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 144.8624 5488.414 4108.696 0.04028533 3.222525 0.5226306
##                     ACF1
## Training set 0.005097859

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. el valor p es la última columna, cuando tiene 1 asterisco es porque es sinificativo con 95% de confianza y *** sería el 100%, así tenga 1 estrella lo vamos a pasar. si todos hubieran salido sin estrella toca descartar el modelo. lo ideal es que todos sean significativo.
coeftest(auto_arima_model_no_seasonal)

## 
## z test of coefficients:
## 
##         Estimate Std. Error z value  Pr(>|z|)    
## ar1    -0.011255   0.126169 -0.0892  0.928919    
## ar2     0.364139   0.115893  3.1420  0.001678 ** 
## ma1    -0.153579   0.089968 -1.7070  0.087816 .  
## ma2    -0.756347   0.087106 -8.6831 < 2.2e-16 ***
## drift 405.625778  71.580859  5.6667 1.456e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

# Mostrar resumen del modelo ajustado
summary(darima_auto)

## Series: train_ts 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.0420  0.3418  -0.0544  -0.6776
## s.e.   0.1317  0.1196   0.0954   0.0848
## 
## sigma^2 = 33338462:  log likelihood = -1530.37
## AIC=3070.74   AICc=3071.15   BIC=3085.86
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 1045.542 5678.818 4410.136 0.7145404 3.431846 0.560974 -0.03008583

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. los graficos de abajo, el de la derecha los datos no están centrados. en este caso el modelo es aceptable pero puede tener errores. El de la izq dice si los residuos están dentro o no dentro de los intervalos, hay 4 líneas o rezagos, no es tan grave, sería grave si todas las líneas están por fuera. esta gráfica se llama correlograma.

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)
## Q* = 69.889, df = 20, p-value = 1.899e-07
## 
## Model df: 4.   Total lags used: 24

Gráfico de Correlograma con algunos residuos por fuera de coeficientes. Prueba aprobada ya que son pocos registros por fuera. Gráfico de distribución normal con cola hacia la izquierda bien distribuido.

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. el test ya lo tomamos como variable previamnete y corresponde a Oct, Nov y Dic, significa que está dentro de muestra.

# 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. LA tendencia lo captura bn el gráfico, centrándose en el punto de quiebre y se encuenra cercano a la línea. se priorizan las últimas 2 observaciones.cuando esté por encima el pronóstico es un sobreajuste y por debajo subajuste

Interpretación modelo automatico (2,1,2):El modelo automático (2,1,2) no parece pronosticar de forma adecuada dentro de prueba ya que no se capturan correctamente los puntos de quiebre. No es un modelo adecuado para predecir futuro.

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

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

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

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

# Mostrar la tabla
print(forecast_table_auto)

##     Tiempo Observado Pronosticado
## 1 2024.750  151683.9     155700.1
## 2 2024.833  158293.7     157380.5
## 3 2024.917  156027.5     157477.3

Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, Enero 2025

Es decir, le sumamos al periodo de prueba una observación más. Se están 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) + 3)

# 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   155700.1
## 2 2024.833   157380.5
## 3 2024.917   157477.3
## 4 2025.000   158350.5
## 5 2025.083   158638.4
## 6 2025.167   159215.7

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

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

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

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

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

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

Identificación automática del modelo SARIMA

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

## Series: train_ts 
## ARIMA(2,1,2)(0,0,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2      ma1      ma2    sma1     drift
##       0.0776  0.2870  -0.1814  -0.7484  0.4582  410.2079
## s.e.  0.1455  0.1362   0.1098   0.1107  0.0855   77.2932
## 
## sigma^2 = 26303091:  log likelihood = -1512.98
## AIC=3039.96   AICc=3040.74   BIC=3061.13

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(2, 1, 2),  # (p,d,q) -> (2,1,2)
                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(2,1,2)(0,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1      ma2    sma1
##       0.0472  0.2504  -0.0957  -0.6712  0.4657
## s.e.  0.1517  0.1405   0.1166   0.1092  0.0840
## 
## sigma^2 = 27592965:  log likelihood = -1516.94
## AIC=3045.88   AICc=3046.46   BIC=3064.03
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 820.2403 5148.873 3854.356 0.5551843 3.010244 0.4902783 -0.0145149

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 al modelo anterior. Al comparar los valores reales VS pronosticados se determina una poca coincidencia. Sigue funcionando mejor el modelo automatico (2,1,2)

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

Análisis 1. Comportamiento de los residuos (Gráfico superior) La serie de residuos oscila alrededor de cero, sin una tendencia clara, lo cual indica que el modelo captura bien la estructura principal de la serie. Sin embargo, hay picos atípicos, especialmente hacia 2020 (coincidiendo con la pandemia), lo que sugiere que el modelo no explica completamente los choques extremos.

Autocorrelación de los residuos (ACF, gráfico inferior izquierdo) La mayoría de los rezagos están dentro de las bandas de confianza (líneas azules), aunque hay algunos picos aislados. Esto indica que los residuos son casi ruido blanco, pero no perfectamente; persisten pequeñas correlaciones que podrían mejorarse con ajustes adicionales.
Distribución de residuos (gráfico inferior derecho) La distribución es aproximadamente normal, centrada en cero, aunque con colas algo más gruesas (kurtosis alta), lo que indica presencia de algunos valores extremos. Esto es aceptable para pronóstico, pero advierte que los intervalos de confianza pueden subestimar el riesgo en eventos atípicos.

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. aunque tenga mejor calificación en los residuos, la siguiente gráfica no se ajusta bn en la caida

Pronóstico del modelo automático SARIMA en el set de prueba-Tabla Predice de mejor manera que elmodelo ARIMA,

El gráfico compara los valores pronosticados (SARIMA) frente a los observados reales para la producción de pollo en Colombia hacia el final de 2024:

Relación entre pronóstico y observado Ambos siguen la misma tendencia ascendente seguida de una leve caída, lo que indica que el modelo capturó bien la dirección general del comportamiento. Sin embargo, el pronóstico subestima los valores reales en el punto máximo (alrededor del periodo 2024.85), mostrando una diferencia aproximada de 2.000 a 2.500 unidades.
Precisión del modelo Aunque hay diferencia en magnitud, el patrón es consistente, lo que indica que el modelo es útil para predicción de tendencias, pero menos preciso en picos de demanda. Esto puede deberse a: Alta volatilidad en los últimos meses de 2024. Posibles factores externos no modelados (ej. costos de insumos, políticas agrícolas, importaciones de cereales).

# 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  151683.9     153599.7
## 2 2024.833  158293.7     158675.9
## 3 2024.917  156027.5     157874.6

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

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

# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente mes (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   153599.7
## 2 2024.833   158675.9
## 3 2024.917   157874.6
## 4 2025.000   154761.1

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

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

Conclusión:

El modelo SARIMA (2,1,2)(0,0,1) es más robusto que el ARIMA simple porque aptura estacionalidad, pero tiende a suavizar quiebres bruscos, lo que limita su precisión en periodos de crisis.

Los residuos confirman que el modelo es estadísticamente válido (sin autocorrelación), pero en la práctica el sector está sujeto a choques no modelados (políticos, sanitarios, climáticos.

El pronóstico sugiere que el sector de pollo se mantendría estable a inicios de 2025, sin grandes crecimientos, a menos que haya un cambio estructural en costos o demanda.

By Adriana Ruiz / Juan Alvear

Series de tiempo: Extracción de señales y Modelo ARIMA - SARIMA sector pecuario Colombia

Adriana Ruiz / Juan Alvear

2025-08-23

Extracción de señales

Serie ajustada por estacionalidad medianamente suavizada y ajustada a original

Serie ajustada por estacionalidad medianamente suavizada y ajustada a original

Serie temporal de la producción de pollo (Variable 1) frente a su tendencia obtenida por una descomposición estacional. En los último años la tendencia es estable, con una ligera caida finalizando el año 2024.

La caida reflejada entre el año 2022 y 2023 debido a la volatilidad de los precios de los productos básicos, el aumento de la inflación, el fortalecimiento del dólar y el contexto de la guerra en Ucrania que afectó los precios y la oferta de alimentos a nivel mundial.

Modelo ARIMA

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

Gráfico de Correlograma con algunos residuos por fuera de coeficientes. Prueba aprobada ya que son pocos registros por fuera. Gráfico de distribución normal con cola hacia la izquierda bien distribuido.

Modelo SARIMA automático