avance1
Juan Jose Armenta Mindiola
2025-11-04
Mercado de productos lácteos y variación del volúmen
Los consumidores colombianos priorizan cada vez más la salud y la nutrición, impulsando la demanda de productos lácteos ricos en proteínas, calcio, vitaminas y probióticos. Los productos, como la leche fortificada, los yogures ricos en proteínas y las bebidas lácteas funcionales, están ganando popularidad a medida que los padres y los entusiastas del fitness buscan opciones ricas en nutrientes.
Mercado de Productos Lácteos en Colombia – Por Tipo de Producto (Leche, Queso, Postre, Yogur, Mantequilla, Otros); Por Canal de Distribución (Hipermercados y Supermercados, Tiendas Especializadas, Venta minorista en Línea, Tiendas de Conveniencia y Otros); Dinámica del Mercado (2025-2034) y Panorama Competitivo
El valor de la industria de productos lácteos en Colombia está creciendo con la afluencia de iniciativas del gobierno colombiano para promover una mejor calidad de la leche a través de la selección genética, el manejo de pastos, el pastoreo rotativo y las mejoras nutricionales. Por ejemplo, el Plan Nacional de Desarrollo Lechero es la estrategia nacional de Colombia para impulsar la productividad, la calidad y la sostenibilidad del sector lácteo a través del mejoramiento genético, un mejor manejo de los pastos, la capacitación de los agricultores y una mayor infraestructura.
Tendencias e impulsores
Mecanización y adopción de tecnología agrícola moderna: El mercado de productos lácteos en Colombia está ganando terreno a través de la modernización agrícola, ya que incluye el uso de sistemas de ordeño automatizados, software de gestión de rebaños y logística de cadena de frío para aumentar la eficiencia y reducir el desperdicio.
Aumento de los ingresos disponibles y premiumización: La creciente clase media de Colombia y el aumento de los niveles de ingresos permiten a los consumidores cambiar a ofertas de productos lácteos premium, como quesos artesanales, yogures especiales y productos enriquecidos.
Innovación en alternativas lácteas y leches vegetales: Las alternativas lácteas de origen vegetal, como las leches de almendras, soja, avena y coco, junto con formulaciones locales que utilizan plátano o coco, están ampliando la participación en el mercado de productos lácteos en Colombia.
Urbanización y cambios en el estilo de vida: La aceleración de la urbanización altera los hábitos alimenticios a medida que más colombianos adoptan estilos de vida ocupados y orientados a la conveniencia. Los envases de Onthego, como tetrapacks, latas de aluminio y bocadillos lácteos listos para comer, abordan las necesidades de los consumidores urbanos.
Interés en comidas listas para comer: La creciente prominencia de las comidas listas para comer convenientes y nutritivas está impulsando significativamente el mercado de productos lácteos en Colombia.
Descripción de la base de datos:
La base de datos contiene 255 observaciones y 11 variables que registra la comercialización del producto de Leche en polvo descremada en dos departamentos de Colombia, con información mensual desde 2007 hasta septiembre del 2018 tomada de Agronet. Cada fila corresponde a una observación regional e incluye las variables: Departamento, Producto, Fecha (mes y año), mes, año, Precio (promedio), Unidad Precio, Volúmen (cantidad comercializada), Unidad Volúmen , Variación Precio y Variación Volúmen . La estructura de la tabla permite realizar análisis de tendencias temporales, comparaciones interdepartamentales y estudios sobre la dinámica precio–cantidad en el sector lácteo colombiano, se trabajará con el departamento de Antioquia.
comercioplacteos<- read_excel("Datos305 (1).xlsx")
comercioplacteos$Año <- as.numeric(format(as.Date(comercioplacteos$Año), "%Y"))#dimensiones
dimensiones <- dim(comercioplacteos)
print(dimensiones)## [1] 255 11#Verificar nombres de las variables
variablesnombres <- names(comercioplacteos)
print(variablesnombres)## [1] "Departamento" "Producto" "Año"
## [4] "Mes" "Fecha" "Precio"
## [7] "Unidad Precio" "Volúmen" "Unidad Volúmen"
## [10] "Variación Precio (%)" "Variación Volúmen (%)"#Tipos de datos
str(comercioplacteos)## tibble [255 × 11] (S3: tbl_df/tbl/data.frame)
## $ Departamento : chr [1:255] "Cundinamarca" "Antioquia" "Cundinamarca" "Antioquia" ...
## $ Producto : chr [1:255] "Leche en polvo descremada" "Leche en polvo descremada" "Leche en polvo descremada" "Leche en polvo descremada" ...
## $ Año : num [1:255] 2007 2007 2007 2007 2007 ...
## $ Mes : POSIXct[1:255], format: "2007-01-01" "2007-02-01" ...
## $ Fecha : chr [1:255] "Enero de 2007" "Febrero de 2007" "Febrero de 2007" "Marzo de 2007" ...
## $ Precio : num [1:255] 7712 8375 10078 9246 11637 ...
## $ Unidad Precio : chr [1:255] "$/Kg" "$/Kg" "$/Kg" "$/Kg" ...
## $ Volúmen : num [1:255] 28963 291313 21315 233038 100745 ...
## $ Unidad Volúmen : chr [1:255] "Kg" "Kg" "Kg" "Kg" ...
## $ Variación Precio (%) : num [1:255] -32.6 8.6 -13.4 10.4 50.9 ...
## $ Variación Volúmen (%): num [1:255] 90.3 905.8 -78.8 -20 247.8 ...Descriptivo a la base de datos por el departamento de antioquia
# Filtrar solo Antioquia
antioquia_data <- comercioplacteos %>% filter(Departamento == "Antioquia")
describe(antioquia_data)## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf## vars n mean sd median trimmed mad
## Departamento* 1 127 1.00 0.00 1.00 1.00 0.00
## Producto* 2 127 1.00 0.00 1.00 1.00 0.00
## Año 3 127 2012.06 3.34 2012.00 2011.97 4.45
## Mes 4 127 NaN NA NA NaN NA
## Fecha* 5 127 64.00 36.81 64.00 64.00 47.44
## Precio 6 127 13447.45 2053.87 13337.00 13404.09 1982.24
## Unidad Precio* 7 127 1.00 0.00 1.00 1.00 0.00
## Volúmen 8 127 10571.12 34360.88 2953.38 3763.22 1973.04
## Unidad Volúmen* 9 127 1.00 0.00 1.00 1.00 0.00
## Variación Precio (%) 10 127 1.16 10.18 -0.02 0.30 6.24
## Variación Volúmen (%) 11 127 161.78 838.64 0.00 6.61 51.90
## min max range skew kurtosis se
## Departamento* 1.00 1.00 0.00 NaN NaN 0.00
## Producto* 1.00 1.00 0.00 NaN NaN 0.00
## Año 2007.00 2018.00 11.00 0.20 -1.07 0.30
## Mes Inf -Inf -Inf NA NA NA
## Fecha* 1.00 127.00 126.00 0.00 -1.23 3.27
## Precio 8375.00 18002.40 9627.40 0.19 -0.35 182.25
## Unidad Precio* 1.00 1.00 0.00 NaN NaN 0.00
## Volúmen 1.00 291313.00 291312.00 6.54 46.12 3049.04
## Unidad Volúmen* 1.00 1.00 0.00 NaN NaN 0.00
## Variación Precio (%) -26.80 38.75 65.55 0.99 2.54 0.90
## Variación Volúmen (%) -98.51 7500.00 7598.51 6.99 52.84 74.42Por las variables de estudio
Histograma para volumen
options(scipen=999)
# Histograma con curva de densidad
ggplot(data.frame(Volúmen = antioquia_data$Volúmen), aes(x = Volúmen)) +
geom_histogram(aes(y = ..density..), bins = 30, fill = "orange", color = "black") +
geom_density(color = "red") +
labs(title = "Histograma con curva de densidad - Volúmen", x = "Volúmen", y = "Kg")## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
datos_volumen <- antioquia_data$Volúmen
# Ajustar las distribuciones
fit_normal <- fitdistr(datos_volumen, "normal")
fit_lognormal <- fitdistr(datos_volumen, "lognormal")
fit_exp <- fitdistr(datos_volumen, "exponential")
# Mostrar AICs y BICs
cat("=== COMPARACIÓN DE AIC Y BIC ===\n")## === COMPARACIÓN DE AIC Y BIC ===cat("Normal - AIC:", AIC(fit_normal), "BIC:", BIC(fit_normal), "\n")## Normal - AIC: 3016.354 BIC: 3022.042cat("Log-Normal - AIC:", AIC(fit_lognormal), "BIC:", BIC(fit_lognormal), "\n")## Log-Normal - AIC: 2503.451 BIC: 2509.139cat("Exponencial - AIC:", AIC(fit_exp), "BIC:", BIC(fit_exp), "\n")## Exponencial - AIC: 2609.534 BIC: 2612.378# Identificar la mejor distribución por AIC y BIC
aic_values <- c(
Normal = AIC(fit_normal),
LogNormal = AIC(fit_lognormal),
Exponencial = AIC(fit_exp)
)
bic_values <- c(
Normal = BIC(fit_normal),
LogNormal = BIC(fit_lognormal),
Exponencial = BIC(fit_exp)
)
mejor_dist_aic <- names(which.min(aic_values))
mejor_dist_bic <- names(which.min(bic_values))
cat("\n=== RESULTADOS ===\n")##
## === RESULTADOS ===cat("MEJOR DISTRIBUCIÓN por AIC:", mejor_dist_aic, "con AIC =", round(min(aic_values), 2), "\n")## MEJOR DISTRIBUCIÓN por AIC: LogNormal con AIC = 2503.45cat("MEJOR DISTRIBUCIÓN por BIC:", mejor_dist_bic, "con BIC =", round(min(bic_values), 2), "\n")## MEJOR DISTRIBUCIÓN por BIC: LogNormal con BIC = 2509.14# Crear secuencia de x para las curvas
x_seq <- seq(min(datos_volumen), max(datos_volumen), length.out = 1000)
# Data frame para las curvas de densidad
curvas_df <- data.frame(
x = rep(x_seq, 3),
y = c(
dnorm(x_seq, mean = fit_normal$estimate["mean"], sd = fit_normal$estimate["sd"]),
dlnorm(x_seq, meanlog = fit_lognormal$estimate["meanlog"], sdlog = fit_lognormal$estimate["sdlog"]),
dexp(x_seq, rate = fit_exp$estimate["rate"])
),
Distribucion = rep(c("Normal", "Log-Normal", "Exponencial"), each = length(x_seq))
)
# Crear tabla de resultados para mostrar en el gráfico
resultados_texto <- paste(
"Mejor por AIC: ", mejor_dist_aic, " (AIC = ", round(min(aic_values), 2), ")\n",
"Mejor por BIC: ", mejor_dist_bic, " (BIC = ", round(min(bic_values), 2), ")",
sep = ""
)
ggplot() +
# Histograma de datos observados
geom_histogram(data = data.frame(Volúmen = datos_volumen),
aes(x = Volúmen, y = ..density..),
bins = 30, fill = "lightblue", color = "black", alpha = 0.6) +
# Densidad de datos observados
geom_density(data = data.frame(Volúmen = datos_volumen),
aes(x = Volúmen), color = "black", linewidth = 1.5, linetype = "solid") +
# Curvas de distribuciones teóricas
geom_line(data = curvas_df,
aes(x = x, y = y, color = Distribucion),
linewidth = 1.2, alpha = 0.8) +
# Anotación con resultados
annotate("text", x = Inf, y = Inf,
label = resultados_texto,
hjust = 1.1, vjust = 1.1,
size = 4, color = "darkred", fontface = "bold") +
labs(title = "Comparación de Ajustes de Distribuciones - Volúmen",
subtitle = "Evaluación por AIC y BIC",
x = "Volúmen (Kg)",
y = "Densidad") +
scale_color_manual(values = c("Normal" = "#ACE0B4",
"Log-Normal" = "lightblue",
"Exponencial" = "#E6B9D9")) +
theme_minimal() +
theme(legend.position = "bottom",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11))
La distribución del volumen distribuido es muy asimétrica, esto revela que existen meses con ventas más altas que otras. El ajuste más adecuado a una distribución log-normal confirma que los datos tienen un sesgo positivo demostrando que el volumen se ven afectados por otros elementos.
Histograma para precio
densidad_precio <- density(antioquia_data$Precio, na.rm = TRUE)
hist_dens_precio <- ggplot(comercioplacteos, aes(x = Precio)) +
geom_histogram(aes(y = ..density..),
bins = 30,
fill = "skyblue",
color = "black",
alpha = 0.7) +
geom_density(color = "red", linewidth = 1) +
geom_vline(aes(xintercept = mean(Precio, na.rm = TRUE)),
color = "blue",
linewidth = 1,
linetype = "dashed") +
geom_vline(aes(xintercept = median(Precio, na.rm = TRUE)),
color = "green",
linewidth = 1,
linetype = "dashed") +
labs(title = "Histograma de Precios con Densidad",
x = "Precio ($/Kg)",
y = "Densidad") +
theme_minimal() +
annotate("text",
x = max(comercioplacteos$Precio, na.rm = TRUE) * 0.7,
y = max(densidad_precio$y) * 0.8,
label = paste("Media:", round(mean(comercioplacteos$Precio, na.rm = TRUE), 2)),
color = "blue") +
annotate("text",
x = max(comercioplacteos$Precio, na.rm = TRUE) * 0.7,
y = max(densidad_precio$y) * 0.7,
label = paste("Mediana:", round(median(comercioplacteos$Precio, na.rm = TRUE), 2)),
color = "green")
print(hist_dens_precio)
datos_precio <- comercioplacteos %>%
mutate(Precio = as.numeric(gsub("[\\$,]", "", Precio))) %>%
pull(Precio) %>%
na.omit()
cat("\n=== AJUSTE DE DISTRIBUCIONES PARA PRECIO ===\n")##
## === AJUSTE DE DISTRIBUCIONES PARA PRECIO ===# Lista para almacenar los ajustes
distribuciones <- list()
# Distribución Normal
tryCatch({
fit_normal <- fitdistr(datos_precio, "normal")
cat("Normal - AIC:", AIC(fit_normal), "BIC:", BIC(fit_normal), "\n")
distribuciones[["Normal"]] <- fit_normal
}, error = function(e) {
cat("Normal: No se pudo ajustar\n")
})## Normal - AIC: 4615.427 BIC: 4622.509# Distribución Log-Normal
tryCatch({
fit_lognormal <- fitdistr(datos_precio, "lognormal")
cat("Log-Normal - AIC:", AIC(fit_lognormal), "BIC:", BIC(fit_lognormal), "\n")
distribuciones[["LogNormal"]] <- fit_lognormal
}, error = function(e) {
cat("Log-Normal: No se pudo ajustar\n")
})## Log-Normal - AIC: 4620.839 BIC: 4627.922# Distribución Gamma
tryCatch({
fit_gamma <- fitdistr(datos_precio, "gamma")
cat("Gamma - AIC:", AIC(fit_gamma), "BIC:", BIC(fit_gamma), "\n")
distribuciones[["Gamma"]] <- fit_gamma
}, error = function(e) {
cat("Gamma: No se pudo ajustar\n")
})## Gamma - AIC: 4617.237 BIC: 4624.32# Distribución Weibull
tryCatch({
fit_weibull <- fitdistr(datos_precio, "weibull")
cat("Weibull - AIC:", AIC(fit_weibull), "BIC:", BIC(fit_weibull), "\n")
distribuciones[["Weibull"]] <- fit_weibull
}, error = function(e) {
cat("Weibull: No se pudo ajustar\n")
})## Weibull - AIC: 4625.155 BIC: 4632.237# Recolectar AICs y BICs
aic_values <- c()
bic_values <- c()
for(dist_name in names(distribuciones)) {
aic_values[dist_name] <- AIC(distribuciones[[dist_name]])
bic_values[dist_name] <- BIC(distribuciones[[dist_name]])
}
# Identificar las mejores distribuciones
if(length(aic_values) > 0) {
mejor_dist_aic <- names(which.min(aic_values))
mejor_dist_bic <- names(which.min(bic_values))
cat("\n=== RESULTADOS ===\n")
cat("MEJOR DISTRIBUCIÓN por AIC:", mejor_dist_aic, "con AIC =", round(min(aic_values), 2), "\n")
cat("MEJOR DISTRIBUCIÓN por BIC:", mejor_dist_bic, "con BIC =", round(min(bic_values), 2), "\n")
} else {
stop("No se pudo ajustar ninguna distribución a los datos de precio")
}##
## === RESULTADOS ===
## MEJOR DISTRIBUCIÓN por AIC: Normal con AIC = 4615.43
## MEJOR DISTRIBUCIÓN por BIC: Normal con BIC = 4622.51# Crear secuencia de x para las curvas
x_seq <- seq(min(datos_precio), max(datos_precio), length.out = 1000)
# Data frame para las curvas de densidad
curvas_df <- data.frame()
# Colores pastel para las distribuciones
colores_pastel <- c(
"Normal" = "#FFB6C1", # Rosa pastel
"LogNormal" = "#87CEFA", # Azul pastel
"Gamma" = "#98FB98", # Verde pastel
"Weibull" = "#DDA0DD" # Violeta pastel
)
# Generar curvas para cada distribución ajustada
if("Normal" %in% names(distribuciones)) {
curvas_df <- rbind(curvas_df, data.frame(
x = x_seq,
y = dnorm(x_seq, mean = distribuciones[["Normal"]]$estimate["mean"],
sd = distribuciones[["Normal"]]$estimate["sd"]),
Distribucion = "Normal"
))
}
if("LogNormal" %in% names(distribuciones)) {
curvas_df <- rbind(curvas_df, data.frame(
x = x_seq,
y = dlnorm(x_seq, meanlog = distribuciones[["LogNormal"]]$estimate["meanlog"],
sdlog = distribuciones[["LogNormal"]]$estimate["sdlog"]),
Distribucion = "LogNormal"
))
}
if("Gamma" %in% names(distribuciones)) {
curvas_df <- rbind(curvas_df, data.frame(
x = x_seq,
y = dgamma(x_seq, shape = distribuciones[["Gamma"]]$estimate["shape"],
rate = distribuciones[["Gamma"]]$estimate["rate"]),
Distribucion = "Gamma"
))
}
if("Weibull" %in% names(distribuciones)) {
curvas_df <- rbind(curvas_df, data.frame(
x = x_seq,
y = dweibull(x_seq, shape = distribuciones[["Weibull"]]$estimate["shape"],
scale = distribuciones[["Weibull"]]$estimate["scale"]),
Distribucion = "Weibull"
))
}
# Texto para la anotación
resultados_texto <- paste(
"Mejor por AIC: ", mejor_dist_aic, " (AIC = ", round(aic_values[mejor_dist_aic], 2), ")\n",
"Mejor por BIC: ", mejor_dist_bic, " (BIC = ", round(bic_values[mejor_dist_bic], 2), ")",
sep = ""
)
# Crear el histograma con fondo blanco
ggplot() +
# Histograma con color pastel
geom_histogram(data = data.frame(Precio = datos_precio),
aes(x = Precio, y = ..density..),
bins = 30, fill = "#E6E6FA", color = "#D8BFD8", alpha = 0.7) +
# Densidad de datos observados
geom_density(data = data.frame(Precio = datos_precio),
aes(x = Precio), color = "#2F4F4F", linewidth = 1.5) +
# Curvas de distribuciones teóricas con colores pastel
geom_line(data = curvas_df,
aes(x = x, y = y, color = Distribucion),
linewidth = 1.2, alpha = 0.8) +
# Anotación con resultados
annotate("text", x = Inf, y = Inf,
label = resultados_texto,
hjust = 1.1, vjust = 1.1,
size = 4, color = "darkblue", fontface = "bold") +
labs(title = "Ajuste de Distribuciones para Precio",
subtitle = "Comparación por AIC y BIC",
x = "Precio ($/Kg)",
y = "Densidad",
color = "Distribución") +
scale_color_manual(values = colores_pastel) +
theme_minimal() +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11),
panel.background = element_rect(fill = "white", color = NA),
plot.background = element_rect(fill = "white", color = NA)
)
Los precios se adecuan a una distribución normal, con media y mediana cercanas entre ellas dando evidencia de un comportamiento estable y simétrico a lo largo del tiempo. Esto indica que la leche en polvo descremada en Antioquia está influido por el mercado.
Boxplot para volúmen
# Calcular media y mediana para Antioquia
media_volumen_antioquia <- mean(antioquia_data$Volúmen, na.rm = TRUE)
mediana_volumen_antioquia <- median(antioquia_data$Volúmen, na.rm = TRUE)
# Boxplot de Volumen para Antioquia
boxplot_volumen_antioquia <- ggplot(antioquia_data, aes(x = "", y = Volúmen)) +
# Puntos de dispersión
geom_jitter(width = 0.2, alpha = 0.6, color = "#4682B4", size = 1.5) +
# Boxplot
geom_boxplot(fill = "#FFA500", color = "black", width = 0.3,
outlier.color = "red", outlier.shape = 16, alpha = 0.7) +
# Media
stat_summary(fun = mean, geom = "point", shape = 23, size = 4,
fill = "blue", aes(x = "")) +
# Mediana
geom_hline(yintercept = mediana_volumen_antioquia, linetype = "dashed",
color = "green", linewidth = 1) +
# Escala logarítmica
scale_y_log10() +
theme_minimal(base_size = 14) +
labs(title = "Distribución del Volumen - Antioquia",
subtitle = paste("Media (rombo azul):", round(media_volumen_antioquia, 2),
"| Mediana (línea verde):", round(mediana_volumen_antioquia, 2)),
y = "Volumen (log10)",
x = "") +
coord_flip() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
annotate("text", x = 1.2, y = media_volumen_antioquia,
label = "Media", color = "blue", size = 4) +
annotate("text", x = 1.2, y = mediana_volumen_antioquia,
label = "Mediana", color = "green", size = 4)
print(boxplot_volumen_antioquia)
El boxplot indica que existe alta asimetría en la distribución con una mediana muy cercana y una media alta. La presencia de una cantidad considerable de outliers por encima del bigote superior indica meses con volumenes altos, está distribución suguiere una comercialización irregular.
Boxplot para precio
# Calcular media y mediana de precio para Antioquia
media_precio_antioquia <- mean(antioquia_data$Precio, na.rm = TRUE)
mediana_precio_antioquia <- median(antioquia_data$Precio, na.rm = TRUE)
# Boxplot de Precio para Antioquia
boxplot_precio_antioquia <- ggplot(antioquia_data, aes(x = "", y = Precio)) +
# Puntos de dispersión
geom_jitter(width = 0.2, alpha = 0.6, color = "#4682B4", size = 1.5) +
# Boxplot
geom_boxplot(fill = "lightgreen", alpha = 0.7, outlier.color = "red", outlier.size = 2) +
# Media
stat_summary(fun = mean, geom = "point", shape = 23, size = 4, fill = "blue") +
# Línea de mediana
geom_hline(yintercept = mediana_precio_antioquia, linetype = "dashed", color = "darkgreen", linewidth = 1) +
labs(title = "Distribución de Precios - Antioquia",
subtitle = paste("Media (rombo azul):", round(media_precio_antioquia, 2),
"| Mediana (línea verde):", round(mediana_precio_antioquia, 2)),
x = "",
y = "Precio ($/Kg)") +
theme_minimal() +
coord_flip() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 11),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank()
) +
# Anotaciones para media y mediana
annotate("text", x = 1.2, y = media_precio_antioquia,
label = "Media", color = "blue", size = 4, fontface = "bold") +
annotate("text", x = 1.2, y = mediana_precio_antioquia,
label = "Mediana", color = "darkgreen", size = 4, fontface = "bold")
print(boxplot_precio_antioquia)
El boxplot del precio muestra una distribución más equilibrada con la media y la mediana cercanas entre sí, la caja centrada sin valores extremos resaltantes, lo que refleja estabilidad en el comportamiento de los precios a lo largo del tiempo.
scatterplot
antioquia_data <- comercioplacteos %>% filter(Departamento == "Antioquia")
scatter_plot_antioquia <- ggplot(antioquia_data, aes(x = Volúmen, y = Precio)) +
geom_point(alpha = 0.7, color = "purple", size = 2.5) +
geom_smooth(method = "lm", color = "red", se = TRUE, fill = "pink", alpha = 0.3) +
labs(title = "Relación entre Precio y Volúmen - Antioquia",
x = "Volúmen (Kg)",
y = "Precio ($/Kg)") +
theme_minimal()
print(scatter_plot_antioquia)## `geom_smooth()` using formula = 'y ~ x'
Correlación
correlacion <- cor(antioquia_data$Precio, antioquia_data$Volúmen, use = "complete.obs");correlacion## [1] -0.2941325Se observa una correlación negativa moderada entre el volumen y el precio, lo que señala que, en términos generales, los aumentos de los precios se vinculan con la reducción de la cantidad comercializada demostrada por el scatterplot y la correlación.
Series de tiempo
antioquia <- comercioplacteos %>%
filter(Departamento == "Antioquia")
# Filtrar Antioquia y crear secuencia temporal
antioquia <- comercioplacteos %>%
filter(Departamento == "Antioquia") %>%
mutate(
Año = as.numeric(substr(Año, 1, 4)),
Mes_num = as.numeric(substr(Mes, 6, 7)),
Fecha = make_date(year = Año, month = Mes_num, day = 1)
) %>%
arrange(Fecha) %>%
mutate(Indice = 1:n()) # Cr# Crear etiquetas para el eje X mostrando los años
etiquetas_años <- antioquia %>%
group_by(Año) %>%
summarise(Posicion = min(Indice))
ggplot(antioquia, aes(x = Indice, y = Precio)) +
geom_line(color = "steelblue", size = 1, alpha = 0.7) +
geom_point(color = "darkblue", size = 2, alpha = 0.8) +
scale_x_continuous(
breaks = etiquetas_años$Posicion,
labels = etiquetas_años$Año,
expand = expansion(mult = c(0.02, 0.02)) # Reduce espacio en los extremos
) +
scale_y_continuous(
labels = scales::dollar_format(prefix = "$", suffix = "/Kg"),
expand = expansion(mult = c(0.05, 0.1)) # Espacio para etiquetas
) +
labs(
title = "Serie de Tiempo: Precio Mensual - Antioquia",
subtitle = "Datos mensuales del precio por kilogramo",
x = "Año",
y = "Precio"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5, color = "gray40"),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
axis.text.y = element_text(size = 9),
panel.grid.major = element_line(color = "gray90", linewidth = 0.2),
panel.grid.minor = element_blank(),
plot.margin = margin(1, 1, 1, 1, "cm")
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
La serie temporal de precios se observan patrones idénticos durante los años del 2007 hasta el 2014, siendo un comportamiento que refleja la inflación acumulada y posibles aumentos o disminuciones de los costos de producción.
# Filtrar Antioquia y crear secuencia temporal
antioquia_volumen <- comercioplacteos %>%
filter(Departamento == "Antioquia") %>%
mutate(
Año = as.numeric(substr(Año, 1, 4)),
Mes_num = as.numeric(substr(Mes, 6, 7)),
Fecha = make_date(year = Año, month = Mes_num, day = 1)
) %>%
arrange(Fecha) %>%
mutate(Indice = 1:n())
ggplot(antioquia_volumen, aes(x = Fecha, y = Volúmen)) +
geom_line(color = "forestgreen", size = 1, alpha = 0.8) +
geom_point(color = "forestgreen", size = 1, alpha = 0.4) + # Puntos más pequeños y transparentes
scale_x_date(
date_breaks = "1 year",
date_labels = "%Y"
) +
labs(
title = "Serie de Tiempo: Volumen Mensual - Antioquia",
x = "Año",
y = "Volumen (Kg)"
) +
theme_minimal()
La serie de volumen es contraria a la de precio con patrones con cambios altos y picos extremos dispersos a lo largo de la serie, no se puede notar tendencias claras de crecimiento y decrecimiento.
# Filtrar Antioquia y crear series de tiempo
antioquia_ts <- comercioplacteos %>%
filter(Departamento == "Antioquia") %>%
mutate(
Año = as.numeric(substr(Año, 1, 4)),
Mes_num = as.numeric(substr(Mes, 6, 7)),
Fecha = make_date(year = Año, month = Mes_num, day = 1)
) %>%
arrange(Fecha)
# Crear series de tiempo mensuales
precio_ts <- ts(antioquia_ts$Precio, start = c(2007, 1), frequency = 12)
volumen_ts <- ts(antioquia_ts$Volúmen, start = c(2007, 1), frequency = 12)Descomposición de la serie precio
# Descomposición aditiva del precio
descomposicion_precio_aditiva <- decompose(precio_ts, type = "additive")
# Graficar SIN el argumento main que causa conflicto
plot(descomposicion_precio_aditiva)
# Descomposición multiplicativa del precio
descomposicion_precio_multiplicativa <- decompose(precio_ts, type = "multiplicative")
# Graficar SIN el argumento main
plot(descomposicion_precio_multiplicativa)
La descomposición aditiva muestra una tendencia creciente hasta el 2015, con componente estacional regular; la descomposición multiplicativa tiene patrones similares demostrando que la estacionariedad es proporcional a la serie.
Descomposición de la serie volúmen
# Descomposición aditiva del volumen
descomposicion_volumen_aditiva <- decompose(volumen_ts, type = "additive")
# Graficar SIN el argumento main
plot(descomposicion_volumen_aditiva)
# Descomposición multiplicativa del volumen
descomposicion_volumen_multiplicativa <- decompose(volumen_ts, type = "multiplicative")
# Graficar SIN el argumento main
plot(descomposicion_volumen_multiplicativa)
La descomposición aditiva como la multiplicativa para la serie de volumen, evidencia un componente estacional notorio, con una tendencia menor y residuos de mayor magnitud contrastable a la de precio.
Estacionariedad
La hipótesis nula de estacionariedad alrededor de una tendencia se confronta con el test de Kwiatkowski-Phillips-Schmidt-Shin (KPSS) para comprobar la suposición de estacionariedad. Esta prueba adicional posibilita verificar que las series temporales no tienen elementos de tendencia que pongan en peligro la validez de los modelos de pronóstico (López, 2002).
Test ADF
Hipótesis
\(H_o\): La serie tiene una raíz unitaria.
\(H_1\): La serie no tiene raíz unitaria.
Si el p-valor < 0.05 se rechaza \(H_0\)
Test KPSS
Hipótesis
- \(H_0:\) El proceso es estacionario en tendencia.
- \(H_1:\) La serie tiene una raíz unitaria.
Estacionariedad para precio
# 1. Prueba de Dickey-Fuller Aumentada (ADF)
adf_test <- adf.test(precio_ts, alternative = "stationary")
print("Prueba ADF para Precio:")## [1] "Prueba ADF para Precio:"print(adf_test)##
## Augmented Dickey-Fuller Test
##
## data: precio_ts
## Dickey-Fuller = -2.7907, Lag order = 5, p-value = 0.2479
## alternative hypothesis: stationarykpss.test(precio_ts)## Warning in kpss.test(precio_ts): p-value smaller than printed p-value##
## KPSS Test for Level Stationarity
##
## data: precio_ts
## KPSS Level = 1.983, Truncation lag parameter = 4, p-value = 0.01La prueba ADF original no rechazó la hipótesis nula de raíz unitaria del p-valor=0.2479, pero la prueba KPSS sí lo hizo con respecto a la estacionariedad p-valor=0.01, lo que corroboró que la serie de precios no era estacionaria.
precio_ts %>%ndiffs()## [1] 1# Aplicar primera diferenciación
precio_diff <- diff(precio_ts)
# Prueba ADF en la serie diferenciada
adf_diff <- adf.test(precio_diff, alternative = "stationary")## Warning in adf.test(precio_diff, alternative = "stationary"): p-value smaller
## than printed p-valueprint("Prueba ADF para Precio (Primera Diferenciación):")## [1] "Prueba ADF para Precio (Primera Diferenciación):"print(adf_diff)##
## Augmented Dickey-Fuller Test
##
## data: precio_diff
## Dickey-Fuller = -6.8598, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationarykpss_diff <- kpss.test(precio_diff)## Warning in kpss.test(precio_diff): p-value greater than printed p-valueprint(kpss_diff)##
## KPSS Test for Level Stationarity
##
## data: precio_diff
## KPSS Level = 0.046285, Truncation lag parameter = 4, p-value = 0.1Al momento de realizarse una diferencia, tanto el ADF como el KPSS mostraron la existencia de estacionariedad, indicando que la serie necesitaba de esta para ser estacionaria corraborada por la línea ndiffs.
Se realiza diferencias a la serie precio:
plot(precio_diff, main = "Precio leche en polvo descremada ",
ylab = "Precio ($/Kg)", xlab = "Año", col = "steelblue", lwd = 2)
abline(h = mean(precio_diff, na.rm = TRUE), col = "red", lty = 2, lwd = 2)
grid()
precio_diff%>%ndiffs()## [1] 0Estacionariedad para volúmen
volumen_ts %>%ndiffs()## [1] 0## [1] "Prueba ADF para Volúmen:"##
## Augmented Dickey-Fuller Test
##
## data: volumen_ts
## Dickey-Fuller = -4.0764, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary##
## KPSS Test for Level Stationarity
##
## data: volumen_ts
## KPSS Level = 0.34529, Truncation lag parameter = 4, p-value = 0.1El test ADF aplicado al volumen rechaza la hipótesis nula de raíz unitaria,mientras el KPSS no rechaza la hipótesis nula correspondiente que el proceso es estacionario. Esta coincidencia muestra que ka la serie de volumen ya es estacionadia sin ka necesidad de hacer diferencias como se muestra en ka línea con ndiffs.
ACF y PACF precio_ts, precio_diff
# ---------------------------------------------------------
# ACF y PACF para precio_ts
# ---------------------------------------------------------
# Configurar área de gráficos en 1 fila y 2 columnas
par(mfrow = c(1, 2))
# ACF para precio_ts
acf(precio_ts, main = "ACF - Precio_TS",
lag.max = 10, col = "darkblue")
# PACF para precio_ts
pacf(precio_ts, main = "PACF - Precio_TS",
lag.max = 10, col = "forestgreen")
# Restaurar configuración de gráficos
par(mfrow = c(1, 1))La función de autocorrelación (FAC) presenta una disminución lenta y constante, con valores significativos que perduran a través de varios rezagos. Mientras la función de autocorrección parcial (FACP), presenta un corte abrupto después del primer rezago.
# ---------------------------------------------------------
# ACF y PACF para precio_diff (serie diferenciada)
# ---------------------------------------------------------
par(mfrow = c(1, 2))
# ACF para precio_diff
acf(precio_diff, main = "ACF",
lag.max = 10, col = "red")
# PACF para precio_diff
pacf(precio_diff, main = "PACF ",
lag.max = 10, col = "purple")
# Restaurar configuración de gráficos
par(mfrow = c(1, 1))Después de la diferencia el FAC,muestra un decaimiento rápido y unos pocos rezagos. Mientras el FACP presenta un patrón de decaimiento exponencial con varios rezagos, mostrándose mixta es decir, un modelo ARMA.
# ---------------------------------------------------------
# ACF y PACF para volumen_ts
# ---------------------------------------------------------
par(mfrow = c(1, 2))
# ACF para volumen_ts
acf(volumen_ts, main = "ACF volumen_ts",
lag.max = 10, col = "blue")
# PACF para volumen_ts
pacf(volumen_ts, main = "PACF volumen_ts",
lag.max = 10, col = "black")
# Restaurar configuración de gráficos
par(mfrow = c(1, 1))La función de autocorrelación del volumen presenta cifras relevantes en los rezagos iniciales, las cuales decrecen rápidamente después. En la Autocorrelación parcial presenta un patrón mixto con varios rezagos.
Identificación del modelo
- Identificación para el precio y precio_diff
Modelo para precio sin diferenciar
# ---------------------------------------------------------
# Comparación de modelos ARMA para precio_ts (SIN DIFERENCIACIÓN)
# ---------------------------------------------------------
# Lista de modelos a evaluar (d=0 para serie no estacionaria)
modelos <- list(
# Modelos AR
list(nombre = "AR(1)", orden = c(1, 0, 0)),
list(nombre = "AR(2)", orden = c(2, 0, 0)),
list(nombre = "AR(3)", orden = c(3, 0, 0)),
list(nombre = "AR(4)", orden = c(4, 0, 0)),
list(nombre = "AR(5)", orden = c(5, 0, 0)),
# Modelos MA
list(nombre = "MA(1)", orden = c(0, 0, 1)),
list(nombre = "MA(2)", orden = c(0, 0, 2)),
list(nombre = "MA(3)", orden = c(0, 0, 3)),
list(nombre = "MA(4)", orden = c(0, 0, 4)),
list(nombre = "MA(5)", orden = c(0, 0, 5)),
# Modelos ARMA
list(nombre = "ARMA(1,1)", orden = c(1, 0, 1)),
list(nombre = "ARMA(1,2)", orden = c(1, 0, 2)),
list(nombre = "ARMA(2,1)", orden = c(2, 0, 1)),
list(nombre = "ARMA(2,2)", orden = c(2, 0, 2)),
list(nombre = "ARMA(2,3)", orden = c(2, 0, 3)),
list(nombre = "ARMA(3,1)", orden = c(3, 0, 1)),
list(nombre = "ARMA(3,2)", orden = c(3, 0, 2)),
list(nombre = "ARMA(4,1)", orden = c(4, 0, 1)),
list(nombre = "ARMA(4,2)", orden = c(4, 0, 2))
)
# Crear data frame vacío para resultados
res <- data.frame(
Modelo = character(),
AIC = numeric(),
BIC = numeric(),
stringsAsFactors = FALSE
)
# Ajustar cada modelo sobre precio_ts (SIN DIFERENCIACIÓN)
for (m in modelos) {
tryCatch({
fit <- arima(precio_ts, order = m$orden, include.mean = FALSE,
method = "ML")
res <- rbind(res, data.frame(
Modelo = m$nombre,
AIC = AIC(fit),
BIC = BIC(fit)
))
}, error = function(e) {
cat("Error en modelo", m$nombre, ":", e$message, "\n")
})
}## Warning in arima(precio_ts, order = m$orden, include.mean = FALSE, method =
## "ML"): possible convergence problem: optim gave code = 1
## Warning in arima(precio_ts, order = m$orden, include.mean = FALSE, method =
## "ML"): possible convergence problem: optim gave code = 1# Ordenar por AIC
res <- res[order(res$AIC), ]
rownames(res) <- NULL
# Mostrar resultados
cat("\n===== COMPARACIÓN DE MODELOS ARMA PARA PRECIO_TS (SIN DIFERENCIACIÓN) =====\n")##
## ===== COMPARACIÓN DE MODELOS ARMA PARA PRECIO_TS (SIN DIFERENCIACIÓN) =====print(res)## Modelo AIC BIC
## 1 ARMA(1,2) 2156.155 2167.532
## 2 ARMA(2,2) 2157.568 2171.788
## 3 ARMA(1,1) 2158.207 2166.739
## 4 ARMA(3,1) 2158.347 2172.568
## 5 ARMA(2,3) 2159.580 2176.645
## 6 ARMA(3,2) 2159.591 2176.656
## 7 AR(5) 2160.053 2177.119
## 8 AR(3) 2160.300 2171.676
## 9 ARMA(4,2) 2160.813 2180.722
## 10 AR(4) 2161.245 2175.466
## 11 ARMA(4,1) 2163.883 2180.948
## 12 AR(2) 2170.530 2179.062
## 13 AR(1) 2180.579 2186.268
## 14 ARMA(2,1) 2183.735 2195.111
## 15 MA(5) 2365.555 2382.620
## 16 MA(4) 2384.030 2398.251
## 17 MA(3) 2445.384 2456.761
## 18 MA(2) 2537.346 2545.879
## 19 MA(1) 2624.351 2630.039# Identificar mejor modelo
if(nrow(res) > 0) {
cat("\nMejor modelo sin diferenciacion:", res[1, "Modelo"], "\n")
cat("AIC:", res[1, "AIC"], "\n")
# Advertencia sobre no estacionariedad
cat("\nADVERTENCIA: precio_ts es no estacionaria\n")
cat("Los modelos ARMA (d=0) pueden no ser apropiados\n")
cat("Se recomienda usar modelos ARIMA (d>=1) o diferenciar la serie primero\n")
}##
## Mejor modelo sin diferenciacion: ARMA(1,2)
## AIC: 2156.155
##
## ADVERTENCIA: precio_ts es no estacionaria
## Los modelos ARMA (d=0) pueden no ser apropiados
## Se recomienda usar modelos ARIMA (d>=1) o diferenciar la serie primero# Auto ARIMA para serie no diferenciada
auto_precio <- auto.arima(precio_ts,
d = 0, # No diferenciar - YA es estacionaria
D = 0, # No diferenciación estacional
seasonal = FALSE, # Sin estacionalidad
stepwise = FALSE,
approximation = FALSE,
trace = TRUE)##
## ARIMA(0,0,0) with zero mean : 2780.011
## ARIMA(0,0,0) with non-zero mean : 2300.883
## ARIMA(0,0,1) with zero mean : 2624.447
## ARIMA(0,0,1) with non-zero mean : 2230.793
## ARIMA(0,0,2) with zero mean : 2513.062
## ARIMA(0,0,2) with non-zero mean : 2214.55
## ARIMA(0,0,3) with zero mean : Inf
## ARIMA(0,0,3) with non-zero mean : 2197.362
## ARIMA(0,0,4) with zero mean : 2384.526
## ARIMA(0,0,4) with non-zero mean : 2191.764
## ARIMA(0,0,5) with zero mean : 2366.244
## ARIMA(0,0,5) with non-zero mean : 2184.445
## ARIMA(1,0,0) with zero mean : Inf
## ARIMA(1,0,0) with non-zero mean : 2170.636
## ARIMA(1,0,1) with zero mean : Inf
## ARIMA(1,0,1) with non-zero mean : 2156.332
## ARIMA(1,0,2) with zero mean : Inf
## ARIMA(1,0,2) with non-zero mean : 2154.92
## ARIMA(1,0,3) with zero mean : Inf
## ARIMA(1,0,3) with non-zero mean : 2156.557
## ARIMA(1,0,4) with zero mean : Inf
## ARIMA(1,0,4) with non-zero mean : 2158.599
## ARIMA(2,0,0) with zero mean : Inf
## ARIMA(2,0,0) with non-zero mean : 2165.256
## ARIMA(2,0,1) with zero mean : Inf
## ARIMA(2,0,1) with non-zero mean : 2155.834
## ARIMA(2,0,2) with zero mean : Inf
## ARIMA(2,0,2) with non-zero mean : 2156.454
## ARIMA(2,0,3) with zero mean : Inf
## ARIMA(2,0,3) with non-zero mean : 2158.686
## ARIMA(3,0,0) with zero mean : Inf
## ARIMA(3,0,0) with non-zero mean : 2157.432
## ARIMA(3,0,1) with zero mean : Inf
## ARIMA(3,0,1) with non-zero mean : 2157.186
## ARIMA(3,0,2) with zero mean : Inf
## ARIMA(3,0,2) with non-zero mean : 2158.678
## ARIMA(4,0,0) with zero mean : Inf
## ARIMA(4,0,0) with non-zero mean : 2158.984
## ARIMA(4,0,1) with zero mean : Inf
## ARIMA(4,0,1) with non-zero mean : 2158.818
## ARIMA(5,0,0) with zero mean : Inf
## ARIMA(5,0,0) with non-zero mean : 2158.644
##
##
##
## Best model: ARIMA(1,0,2) with non-zero mean# Resumen del modelo
cat("=== MODELO ARMA PARA PRECIO_ts ===\n")## === MODELO ARMA PARA PRECIO_ts ===summary(auto_precio)## Series: precio_ts
## ARIMA(1,0,2) with non-zero mean
##
## Coefficients:
## ar1 ma1 ma2 mean
## 0.9879 -0.4552 -0.2010 13482.685
## s.e. 0.0152 0.0974 0.1023 1878.383
##
## sigma^2 = 1280310: log likelihood = -1072.21
## AIC=2154.42 AICc=2154.92 BIC=2168.65
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 146.0949 1113.546 829.1434 0.4145067 6.45018 0.545398 -0.04344806# Extraer los órdenes del modelo
p <- auto_precio$arma[1] # orden AR
q <- auto_precio$arma[2] # orden MA
cat("Modelo encontrado: ARMA(", p, ",", q, ")\n", sep = "")## Modelo encontrado: ARMA(1,2)ARMA(1,2)
\(X_t = 13482.685 + 0.9879 x_{t-1} + \epsilon_t - 0.4552 \epsilon_{t-1} - 0.2010 \epsilon_{t-2}\)
El modelo que mejor se ajusta por las dos pruebas realizadas es un ARMA(1,2), para la serie de precio sin realizar la debida diferencia. Un coeficiente perteneciente al modelo AR y dos coeficientes por parte de el modelo MA siendo estas: MA(1) y MA(2), sin teneruna diferencia no se puede decir que este es el mejor modelo porque la serie no es estacionaria.
Precio diferenciada
modelos_arima <- list(
list(nombre = "ARIMA(1,1,1)", orden = c(1, 1, 1)),
list(nombre = "ARIMA(1,1,2)", orden = c(1, 1, 2)),
list(nombre = "ARIMA(1,1,3)", orden = c(1, 1, 3)),
list(nombre = "ARIMA(2,1,1)", orden = c(2, 1, 1)),
list(nombre = "ARIMA(2,1,2)", orden = c(2, 1, 2)),
list(nombre = "ARIMA(2,1,3)", orden = c(2, 1, 3)),
list(nombre = "ARIMA(3,1,1)", orden = c(3, 1, 1)),
list(nombre = "ARIMA(3,1,2)", orden = c(3, 1, 2)),
list(nombre = "ARIMA(3,1,3)", orden = c(3, 1, 3)),
list(nombre = "ARIMA(4,1,1)", orden = c(4, 1, 1)),
list(nombre = "ARIMA(4,1,2)", orden = c(4, 1, 2)),
list(nombre = "ARIMA(4,1,3)", orden = c(4, 1, 3)),
list(nombre = "ARIMA(0,1,1)", orden = c(0, 1, 1)),
list(nombre = "ARIMA(0,1,2)", orden = c(0, 1, 2)),
list(nombre = "ARIMA(0,1,3)", orden = c(0, 1, 3)),
list(nombre = "ARIMA(1,1,0)", orden = c(1, 1, 0)),
list(nombre = "ARIMA(2,1,0)", orden = c(2, 1, 0)),
list(nombre = "ARIMA(3,1,0)", orden = c(3, 1, 0)),
list(nombre = "ARIMA(4,1,0)", orden = c(4, 1, 0))
)
# Crear data frame vacío para resultados
res_arima <- data.frame(
Modelo = character(),
AIC = numeric(),
BIC = numeric(),
Coeficientes = character(),
stringsAsFactors = FALSE
)
# Lista para almacenar los modelos ajustados
modelos_ajustados_arima <- list()
# ---------------------------------------------------------
# Ajustar cada modelo ARIMA sobre la serie original (precio_ts)
# ---------------------------------------------------------
for (m in modelos_arima) {
tryCatch({
# Si d=1, incluir drift (para comparar con auto.arima)
fit <- Arima(precio_ts, order = m$orden, include.drift = (m$orden[2] == 1))
# Extraer coeficientes
coef_text <- ""
if (length(fit$coef) > 0) {
coef_names <- names(fit$coef)
coef_values <- round(fit$coef, 4)
coef_text <- paste(paste0(coef_names, "=", coef_values), collapse = ", ")
}
# Guardar modelo y métricas
modelos_ajustados_arima[[m$nombre]] <- fit
res_arima <- rbind(res_arima, data.frame(
Modelo = m$nombre,
AIC = AIC(fit),
BIC = BIC(fit),
Coeficientes = coef_text,
stringsAsFactors = FALSE
))
}, error = function(e) {
cat("Error ajustando modelo", m$nombre, ":", e$message, "\n")
})
}
# ---------------------------------------------------------
# Ordenar por AIC y BIC
# ---------------------------------------------------------
res_arima_AIC <- res_arima[order(res_arima$AIC), ]
res_arima_BIC <- res_arima[order(res_arima$BIC), ]
rownames(res_arima_AIC) <- NULL
rownames(res_arima_BIC) <- NULL
# ---------------------------------------------------------
# TOP 5 MEJORES MODELOS POR AIC
# ---------------------------------------------------------
cat("\n", strrep("=", 60), "\n", sep = "")##
## ============================================================cat("TOP 5 MEJORES MODELOS ARIMA ORDENADOS POR AIC\n")## TOP 5 MEJORES MODELOS ARIMA ORDENADOS POR AICcat(strrep("=", 60), "\n\n", sep = "")## ============================================================for(i in 1:min(5, nrow(res_arima_AIC))) {
cat("POSICIÓN", i, ":", res_arima_AIC$Modelo[i], "\n")
cat("AIC:", round(res_arima_AIC$AIC[i], 3), "| BIC:", round(res_arima_AIC$BIC[i], 3), "\n")
if(res_arima_AIC$Coeficientes[i] != "") {
cat("Coeficientes:", res_arima_AIC$Coeficientes[i], "\n")
}
cat(strrep("-", 60), "\n", sep = "")
}## POSICIÓN 1 : ARIMA(0,1,2)
## AIC: 2131.153 | BIC: 2142.498
## Coeficientes: ma1=-0.4868, ma2=-0.2465, drift=52.6712
## ------------------------------------------------------------
## POSICIÓN 2 : ARIMA(3,1,1)
## AIC: 2131.78 | BIC: 2148.797
## Coeficientes: ar1=0.4592, ar2=0.0091, ar3=0.202, ma1=-1, drift=45.9463
## ------------------------------------------------------------
## POSICIÓN 3 : ARIMA(1,1,1)
## AIC: 2131.842 | BIC: 2143.187
## Coeficientes: ar1=0.3267, ma1=-0.8394, drift=52.196
## ------------------------------------------------------------
## POSICIÓN 4 : ARIMA(1,1,3)
## AIC: 2132.055 | BIC: 2149.073
## Coeficientes: ar1=0.8956, ma1=-1.4239, ma2=0.2396, ma3=0.1843, drift=46.8946
## ------------------------------------------------------------
## POSICIÓN 5 : ARIMA(1,1,2)
## AIC: 2132.404 | BIC: 2146.586
## Coeficientes: ar1=0.8087, ma1=-1.4108, ma2=0.4108, drift=46.1856
## ------------------------------------------------------------# ---------------------------------------------------------
# DETALLE DEL MEJOR MODELO
# ---------------------------------------------------------
cat("\n", strrep("=", 60), "\n", sep = "")##
## ============================================================cat("DETALLE COMPLETO DEL MEJOR MODELO ARIMA:", res_arima_AIC$Modelo[1], "\n")## DETALLE COMPLETO DEL MEJOR MODELO ARIMA: ARIMA(0,1,2)cat(strrep("=", 60), "\n\n", sep = "")## ============================================================mejor_modelo_arima <- modelos_ajustados_arima[[res_arima_AIC$Modelo[1]]]
print(mejor_modelo_arima)## Series: precio_ts
## ARIMA(0,1,2) with drift
##
## Coefficients:
## ma1 ma2 drift
## -0.4868 -0.2465 52.6712
## s.e. 0.0926 0.1023 27.0337
##
## sigma^2 = 1240918: log likelihood = -1061.58
## AIC=2131.15 AICc=2131.48 BIC=2142.5# ---------------------------------------------------------
# SIGNIFICANCIA ESTADÍSTICA DE COEFICIENTES
# ---------------------------------------------------------
cat("\n", strrep("-", 60), "\n", sep = "")##
## ------------------------------------------------------------cat("SIGNIFICANCIA ESTADÍSTICA DE COEFICIENTES\n")## SIGNIFICANCIA ESTADÍSTICA DE COEFICIENTEScat(strrep("-", 60), "\n", sep = "")## ------------------------------------------------------------if(length(mejor_modelo_arima$coef) > 0) {
coef_summary <- data.frame(
Coeficiente = names(mejor_modelo_arima$coef),
Estimación = round(mejor_modelo_arima$coef, 4),
Error_Std = round(sqrt(diag(mejor_modelo_arima$var.coef)), 4),
stringsAsFactors = FALSE
)
coef_summary$t_value <- round(coef_summary$Estimación / coef_summary$Error_Std, 4)
coef_summary$p_value <- round(2 * (1 - pnorm(abs(coef_summary$t_value))), 4)
coef_summary$Significativo <- ifelse(coef_summary$p_value < 0.05, "SÍ", "NO")
print(coef_summary)
} else {
cat("El modelo no tiene coeficientes estimados\n")
}## Coeficiente Estimación Error_Std t_value p_value Significativo
## ma1 ma1 -0.4868 0.0926 -5.2570 0.0000 SÍ
## ma2 ma2 -0.2465 0.1023 -2.4096 0.0160 SÍ
## drift drift 52.6712 27.0337 1.9484 0.0514 NO# ---------------------------------------------------------
# RESUMEN COMPARATIVO GENERAL
# ---------------------------------------------------------
cat("\n", strrep("=", 60), "\n", sep = "")##
## ============================================================cat("RESUMEN COMPARATIVO - MODELOS ARIMA\n")## RESUMEN COMPARATIVO - MODELOS ARIMAcat(strrep("=", 60), "\n", sep = "")## ============================================================cat("Mejor modelo por AIC:", res_arima_AIC$Modelo[1], "(AIC =", round(res_arima_AIC$AIC[1], 2), ")\n")## Mejor modelo por AIC: ARIMA(0,1,2) (AIC = 2131.15 )cat("Mejor modelo por BIC:", res_arima_BIC$Modelo[1], "(BIC =", round(res_arima_BIC$BIC[1], 2), ")\n")## Mejor modelo por BIC: ARIMA(0,1,2) (BIC = 2142.5 )cat("Total de modelos ARIMA evaluados:", nrow(res_arima), "\n")## Total de modelos ARIMA evaluados: 19cat("Rango de AIC:", round(range(res_arima$AIC), 2), "\n")## Rango de AIC: 2131.15 2148.11cat("Rango de BIC:", round(range(res_arima$BIC), 2), "\n")## Rango de BIC: 2142.5 2161.75auto.arima(precio_ts, stepwise = FALSE, approximation = FALSE)## Series: precio_ts
## ARIMA(0,1,2) with drift
##
## Coefficients:
## ma1 ma2 drift
## -0.4868 -0.2465 52.6712
## s.e. 0.0926 0.1023 27.0337
##
## sigma^2 = 1240918: log likelihood = -1061.58
## AIC=2131.15 AICc=2131.48 BIC=2142.5El modelo es ARIMA(0,1,2):
\(X_t = x_{t-1} + \epsilon_t +0.4868 \epsilon_{t-1} +0.2465\epsilon_{t-2}\)
El modelo que mejor se ajusta a la serie de precio con una diferencia es un ARIMA (0,1,2), con un componente integrado correspondiente a la diferencia realizada para que se vuelva estacionaria y dos componentes de medias móviles (MA1, MA2)
- Identificación para el volúmen
# ---------------------------------------------------------
# Comparación de modelos AR, MA, ARMA y ARIMA para volumen_ts
# ---------------------------------------------------------
# Lista de modelos a evaluar
modelos <- list(
# Modelos AR
list(nombre = "AR(1)", orden = c(1, 0, 0)),
list(nombre = "AR(2)", orden = c(2, 0, 0)),
list(nombre = "AR(3)", orden = c(3, 0, 0)),
list(nombre = "AR(4)", orden = c(4, 0, 0)),
list(nombre = "AR(5)", orden = c(5, 0, 0)),
# Modelos MA
list(nombre = "MA(1)", orden = c(0, 0, 1)),
list(nombre = "MA(2)", orden = c(0, 0, 2)),
list(nombre = "MA(3)", orden = c(0, 0, 3)),
list(nombre = "MA(4)", orden = c(0, 0, 4)),
list(nombre = "MA(5)", orden = c(0, 0, 5)),
# Modelos ARMA
list(nombre = "ARMA(1,1)", orden = c(1, 0, 1)),
list(nombre = "ARMA(1,2)", orden = c(1, 0, 2)),
list(nombre = "ARMA(2,1)", orden = c(2, 0, 1)),
list(nombre = "ARMA(2,2)", orden = c(2, 0, 2)),
list(nombre = "ARMA(2,3)", orden = c(2, 0, 3)),
list(nombre = "ARMA(3,1)", orden = c(3, 0, 1)),
list(nombre = "ARMA(3,2)", orden = c(3, 0, 2)),
list(nombre = "ARMA(4,1)", orden = c(4, 0, 1)),
list(nombre = "ARMA(4,2)", orden = c(4, 0, 2))
)
# Crear data frame vacío para resultados
res <- data.frame(
Modelo = character(),
AIC = numeric(),
BIC = numeric(),
stringsAsFactors = FALSE
)
# Ajustar cada modelo sobre precio_diff
for (m in modelos) {
fit <- arima(volumen_ts, order = m$orden, include.mean = FALSE)
res <- rbind(res, data.frame(
Modelo = m$nombre,
AIC = AIC(fit),
BIC = BIC(fit)
))
}
# Ordenar por AIC
res <- res[order(res$AIC), ]
rownames(res) <- NULL
# Mostrar resultados
cat("\n===== Comparación de modelos (ajustados sobre volumen_ts) =====\n")##
## ===== Comparación de modelos (ajustados sobre volumen_ts) =====print(res)## Modelo AIC BIC
## 1 ARMA(2,3) 2846.213 2863.278
## 2 ARMA(4,2) 2847.857 2867.767
## 3 AR(4) 2847.899 2862.120
## 4 AR(5) 2848.898 2865.963
## 5 ARMA(4,1) 2849.238 2866.303
## 6 ARMA(2,1) 2853.037 2864.414
## 7 ARMA(3,2) 2856.916 2873.981
## 8 AR(3) 2857.324 2868.701
## 9 AR(2) 2857.494 2866.026
## 10 ARMA(1,1) 2858.122 2866.654
## 11 ARMA(1,2) 2859.836 2871.213
## 12 AR(1) 2861.207 2866.895
## 13 ARMA(3,1) 2861.388 2875.609
## 14 ARMA(2,2) 2861.988 2876.209
## 15 MA(5) 2873.130 2890.195
## 16 MA(4) 2880.106 2894.327
## 17 MA(3) 2887.814 2899.191
## 18 MA(2) 2899.081 2907.613
## 19 MA(1) 2935.027 2940.716##
## ARIMA(0,0,0) with zero mean : 3025.958
## ARIMA(0,0,0) with non-zero mean : 3016.45
## ARIMA(0,0,1) with zero mean : 2935.124
## ARIMA(0,0,1) with non-zero mean : 2928.286
## ARIMA(0,0,2) with zero mean : 2899.276
## ARIMA(0,0,2) with non-zero mean : 2895.184
## ARIMA(0,0,3) with zero mean : 2888.142
## ARIMA(0,0,3) with non-zero mean : 2885.438
## ARIMA(0,0,4) with zero mean : 2880.602
## ARIMA(0,0,4) with non-zero mean : 2878.395
## ARIMA(0,0,5) with zero mean : 2873.83
## ARIMA(0,0,5) with non-zero mean : 2872.399
## ARIMA(1,0,0) with zero mean : Inf
## ARIMA(1,0,0) with non-zero mean : 2862.379
## ARIMA(1,0,1) with zero mean : 2858.317
## ARIMA(1,0,1) with non-zero mean : 2859.311
## ARIMA(1,0,2) with zero mean : 2860.164
## ARIMA(1,0,2) with non-zero mean : 2861.191
## ARIMA(1,0,3) with zero mean : 2853.949
## ARIMA(1,0,3) with non-zero mean : 2855.065
## ARIMA(1,0,4) with zero mean : 2850.394
## ARIMA(1,0,4) with non-zero mean : 2851.467
## ARIMA(2,0,0) with zero mean : 2857.689
## ARIMA(2,0,0) with non-zero mean : 2858.654
## ARIMA(2,0,1) with zero mean : 2853.365
## ARIMA(2,0,1) with non-zero mean : 2853.931
## ARIMA(2,0,2) with zero mean : 2862.483
## ARIMA(2,0,2) with non-zero mean : 2863.552
## ARIMA(2,0,3) with zero mean : 2846.913
## ARIMA(2,0,3) with non-zero mean : 2847.689
## ARIMA(3,0,0) with zero mean : 2857.652
## ARIMA(3,0,0) with non-zero mean : 2858.563
## ARIMA(3,0,1) with zero mean : 2861.884
## ARIMA(3,0,1) with non-zero mean : 2862.925
## ARIMA(3,0,2) with zero mean : Inf
## ARIMA(3,0,2) with non-zero mean : Inf
## ARIMA(4,0,0) with zero mean : 2848.395
## ARIMA(4,0,0) with non-zero mean : 2849.059
## ARIMA(4,0,1) with zero mean : 2849.938
## ARIMA(4,0,1) with non-zero mean : 2850.555
## ARIMA(5,0,0) with zero mean : 2849.598
## ARIMA(5,0,0) with non-zero mean : 2850.169
##
##
##
## Best model: ARIMA(2,0,3) with zero mean## Series: volumen_ts
## ARIMA(2,0,3) with zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3
## 1.8590 -0.8695 -0.7303 0.2081 0.2823
## s.e. 0.1391 0.1378 0.1738 0.1467 0.1232
##
## sigma^2 = 283489344: log likelihood = -1417.11
## AIC=2846.21 AICc=2846.91 BIC=2863.28
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 481.2057 16502.37 7613.327 -788.9226 1028.041 0.6592819
## ACF1
## Training set -0.01654261ARMA(2,3)
\(X_t = 1.859 x_{t-1} -0.8695 x_{t-2} + \epsilon_t -0.7303 \epsilon_{t-1} + 0.2081 \epsilon_{t-2} + 0.2823 \epsilon_{t-3}\)
El modelo ARMA(2,3) fue el que mejor se ajustó para el volumen. Esto significa que tiene dos coeficiente de un AR y tres de MA.
Referencias
Estadísticas home. (s. f.). https://www.agronet.gov.co/estadistica/Paginas/home.aspx?cod=9
Presno, M.J. y López, A.J. (2002): “Response surface estimates of stationarity tests with a structural break”, Economics Letters, 78, 395-399.
Tróchez, J. & Valencia, M. (2014). Análisis de series temporales en el sector lácteo de Antioquia para detectar efectos de la apertura comercial. Recuperado de: http://hdl.handle.net/20.500.11912/6814.