Series de tiempo: Extracción de señales y Modelo ARIMA para la produccion de Cafe
Natalya Lopez y Pablo Lopez
2025-03-22
INTRODUCCION Fundada en 1959 por la Federación Nacional de Cafeteros de Colombia, Juan Valdez tiene como misión la promoción del café colombiano en los mercados internacionales. La figura de Juan Valdez, representando a un caficultor colombiano, ha sido determinante para humanizar y dar visibilidad al esfuerzo de los productores nacionales, consolidándose con el tiempo como un símbolo mundial del café colombiano y desempeñando un rol fundamental en las campañas publicitarias y estrategias de marketing de la Federación Nacional de Cafeteros.
La marca Juan Valdez® es la única marca de café de relevancia internacional que pertenece a productores de café. En 2002, la Federación Nacional de Cafeteros de Colombia, una organización sin ánimo de lucro que representa a más de 500.000 familias productoras de café, creó Procafecol S.A. con el propósito de explotar la marca Juan Valdez® a través de tiendas de café y negocios de valor agregado. Tras la popularidad y aceptación de la marca en tiendas, se dio inicio a la distribución de los productos Juan Valdez® en otros canales del mercado colombiano, así como a nivel internacional.
Los cafés de Juan Valdez® son sinónimo de alta calidad para los consumidores de café premium. Además, le entregan a los caficultores un mayor valor por la calidad del café e importantes recursos por concepto de regalías provenientes de la venta de cada taza o producto que lleve su firma en el mundo.
El presente documento examina tres variables esenciales para entender cómo la producción y comercialización del café afectan el rendimiento de la empresa Juan Valdez. En él se estudian la producción de café, el precio internacional y las exportaciones, cubriendo el periodo comprendido entre enero de 2012 y diciembre de 2024. Estas variables resultan cruciales para evaluar su repercusión en una de las materias primas más importantes de la compañía, lo que permite detectar tendencias y patrones que fundamenten decisiones estratégicas en áreas de producción, ventas y formulación de estrategias comerciales en un entorno altamente competitivo.
La importancia de las variables analizadas radica en su estrecha vinculación con el desempeño de la empresa, ya que inciden directamente en la producción, distribución y, en última instancia, en su posicionamiento en el mercado global. Este estudio contribuye a comprender de qué manera las variaciones en dichas variables impactan en las operaciones de la organización, posibilitando el diseño de estrategias comerciales y operativas que optimicen su competitividad en el sector cafetero.
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) Cargar base de datos
library(readxl)
Cafe <- read_excel("Cafe.xlsx", col_types = c("date",
"numeric", "numeric", "numeric"))
View(Cafe)PASO INDISPENSABLE: Declarar las variables como series temporales:
Variable 1: Produccion de Cafe en miles de sacos (60kg)
# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(Cafe$variable1, start = c(2012, 1), frequency = 12)
variable1_ts## Jan Feb Mar Apr May Jun Jul
## 2012 535.0000 571.0000 576.0000 580.0000 689.0000 714.0000 668.0000
## 2013 877.0000 625.0000 617.0000 970.0000 937.0000 913.0000 1031.0000
## 2014 1011.0000 874.0000 828.0000 832.0000 1050.0000 944.0000 1236.0000
## 2015 1088.0000 1029.0000 800.0000 924.0000 1165.0000 1240.0000 1463.0000
## 2016 1136.0000 1096.0000 944.0000 1043.0000 1163.0000 1158.0000 1102.0000
## 2017 1275.0000 1293.0000 1020.0000 834.0000 901.0000 1049.0000 1373.0000
## 2018 1131.0000 1212.0000 1037.0000 874.0000 1188.0000 1087.0000 1051.0000
## 2019 1296.0000 1106.0000 914.0000 1031.0000 1115.0000 1211.0000 1317.0000
## 2020 1050.0000 1001.0000 806.0000 744.0000 1186.0000 1362.0000 1310.0000
## 2021 1081.0000 1107.0000 1050.0000 810.0000 616.0000 1052.0000 1209.0000
## 2022 868.0000 928.0000 914.0000 750.0000 1017.0000 951.0000 943.5992
## 2023 868.0000 1025.0000 799.0000 565.8672 806.1955 956.0000 947.0000
## 2024 959.0000 961.0000 865.5326 742.0000 1120.0000 1172.0000 1158.8838
## Aug Sep Oct Nov Dec
## 2012 565.0000 519.0000 653.0000 770.0000 904.0000
## 2013 770.0000 860.0000 1058.0000 1113.0000 1115.0000
## 2014 1151.0000 912.0000 1101.0000 1115.0000 1086.0000
## 2015 1264.0000 1058.0000 1368.0000 1322.0000 1454.0000
## 2016 1189.0000 1034.0000 1395.0000 1653.0000 1319.0000
## 2017 1294.0000 1228.0000 1073.0000 1304.0000 1550.0000
## 2018 1258.0000 1050.0000 1086.0000 1300.0000 1283.0000
## 2019 1119.0000 1088.0000 1369.0000 1506.0000 1680.0000
## 2020 1091.0000 995.0000 1159.0000 1443.0000 1743.0000
## 2021 915.0000 1209.0000 1012.0000 1131.0000 1385.0000
## 2022 949.0000 834.0000 888.0000 1060.0000 981.0000
## 2023 872.0000 849.0000 1157.4585 1282.1066 1220.0000
## 2024 1049.0000 1071.2334 1339.1550 1761.4095 1798.2306Variable 2: Precio externo del cafe en centavos de dolar por libra
# Convertir/declarar el ISE en serie de tiempo mensual
variable2_ts <- ts(Cafe$variable2, start = c(2012, 1), frequency = 12)
variable2_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2012 256.1168 246.3010 226.0661 215.2930 209.5097 185.8843 200.9219 188.5771
## 2013 168.8300 162.5343 161.7474 162.1890 160.0019 148.4153 147.5277 144.2926
## 2014 131.9742 159.6775 210.1481 216.0827 216.4171 195.8677 193.0277 210.4074
## 2015 185.3126 177.0854 155.8761 156.4597 152.3323 151.1303 145.9206 147.1319
## 2016 136.4284 137.1914 144.3697 144.5490 144.5803 153.9430 164.0368 160.7271
## 2017 163.7090 164.8057 159.2329 156.0010 151.2374 146.8607 150.9077 156.4565
## 2018 143.4810 142.3511 139.7903 138.9967 140.5035 139.1723 134.8471 130.9932
## 2019 128.6758 128.9696 125.4671 124.0727 123.1348 132.7700 138.4145 129.8942
## 2020 150.7962 144.1289 157.3050 164.6219 155.8180 149.2377 151.2900 167.6310
## 2021 175.2395 177.2779 179.1478 181.6810 196.2590 206.1064 214.4400 226.9009
## 2022 293.1955 305.5184 287.9070 292.5780 285.1871 301.2423 286.9865 292.5200
## 2023 218.2400 237.0400 228.2448 232.3695 229.0759 215.4367 190.3325 188.7070
## 2024 206.5600 210.1300 207.5470 237.3114 232.3527 249.9400 257.9800 261.0000
## Sep Oct Nov Dec
## 2012 189.2787 183.1594 171.1003 166.0471
## 2013 139.1247 135.1697 125.3610 125.9203
## 2014 208.0950 223.1000 206.7437 193.3329
## 2015 136.4220 142.7784 138.4003 139.6919
## 2016 167.6153 170.5529 179.7767 160.9187
## 2017 151.5747 145.0094 143.8483 142.5123
## 2018 126.5253 138.8045 140.7503 129.6229
## 2019 131.3550 131.8783 143.5175 160.1652
## 2020 169.8400 156.7536 161.6000 169.9900
## 2021 238.3200 257.1257 273.6352 291.9477
## 2022 296.4571 269.4900 225.1000 223.8900
## 2023 186.0000 184.9700 194.4567 206.9235
## 2024 276.2520 280.0396 293.9385 340.5167Variable 3: Exportacion de cafe en miles de dolares
# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable3_ts <- ts(Cafe$variable3, start = c(2012, 1), frequency = 12)
variable3_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2012 197295.8 186693.3 203636.0 121442.5 167643.9 162433.1 151832.1 180244.0
## 2013 177150.5 171569.2 127181.2 152040.1 168982.2 138026.6 159759.0 169039.9
## 2014 142633.8 202636.1 153516.9 224406.6 213123.9 161247.5 229475.0 218584.5
## 2015 288671.5 251209.4 211129.0 182545.4 194717.5 160584.9 239676.4 218630.5
## 2016 163588.4 184225.8 179113.8 281209.6 152586.3 176288.8 128997.8 169424.1
## 2017 232188.4 214883.5 295930.9 217795.8 196637.3 137977.0 230431.8 211174.9
## 2018 255058.1 187220.6 201578.7 205465.6 190580.9 147559.5 166218.3 213792.4
## 2019 260699.1 221412.3 193807.5 186052.9 143552.4 147796.7 205010.7 218999.2
## 2020 249836.3 215605.6 170321.5 157384.0 141864.0 194805.8 270277.4 201786.8
## 2021 255455.2 267999.8 278224.9 243363.4 124318.7 115655.4 309123.1 322175.0
## 2022 332201.4 383477.5 412235.8 327940.2 275625.7 369007.8 394622.5 323525.4
## 2023 250131.8 261153.3 292938.5 212748.0 251607.5 225461.4 227622.0 226253.7
## 2024 231800.6 258106.5 238506.5 261112.5 268005.6 285915.9 291114.4 318115.4
## Sep Oct Nov Dec
## 2012 123549.6 159793.2 144269.2 166159.5
## 2013 152688.2 152603.0 173174.7 191779.4
## 2014 216896.0 252496.3 267402.2 244052.0
## 2015 250151.3 179340.6 191276.3 217993.9
## 2016 196074.9 193648.9 219217.4 429058.8
## 2017 240608.9 215414.3 198250.7 203778.0
## 2018 179237.1 194640.1 191052.0 216342.7
## 2019 167456.2 195326.7 193952.8 242374.7
## 2020 217897.4 173371.9 227635.8 316675.7
## 2021 268316.6 327953.2 316349.7 376919.2
## 2022 339103.7 323641.1 293141.5 348771.0
## 2023 234930.5 188621.9 251091.2 306637.2
## 2024 284093.6 338376.9 325234.9 462940.1Extracción de Señales en Series de Tiempo para la Industria Cafetera - Caso Juan Valdez
Con el propósito de identificar patrones de comportamiento en el mercado que puedan impactar el rendimiento de la empresa Juan Valdez y detectar potenciales oportunidades estratégicas, se ha llevado a cabo un análisis exhaustivo de tres variables fundamentales que influyen de manera significativa en la disponibilidad de una de las materias primas más críticas para la empresa Juan Valdez: el café. Las variables en estudio son:
1)Producción de café: expresada en miles de sacos de 60 kg. 2)Precio del café: medido en centavos de dólar por libra. 3)Exportaciones de café: expresadas en miles de dolares.
La aplicación de técnicas de análisis de series de tiempo permite descomponer estas variables en sus componentes esenciales: tendencia, estacionalidad y ruido aleatorio. Dado que los datos observados reflejan una interaccion de múltiples influencias, el empleo del software estadístico R permite obtener una comprensión detallada y precisa de los factores subyacentes que determinan las fluctuaciones en cada una de las variables analizadas.
El análisis detallado de las señales extraídas y sus interacciones posibilita la identificación de cambios estructurales en la dinámica del mercado, así como la evaluación de su repercusión en la operatividad y la rentabilidad de Juan Valdez. Los hallazgos obtenidos en este estudio constituirán el fundamento para el desarrollo de estrategias orientadas a robustecer la competitividad de la empresa y mitigar los riesgos asociados a la volatilidad del sector.
La fase de extracción de señales se configura como un prerrequisito esencial para la aplicación de modelos de pronóstico, tales como ARIMA y SARIMA, dentro del marco metodológico Box-Jenkins. Estos modelos alcanzan niveles superiores de precisión cuando se incorporan correctamente los componentes subyacentes de las series temporales, especialmente la estacionalidad y la tendencia, lo que permite focalizar el poder predictivo en las variaciones sistemáticas y reducir la incidencia de errores derivados de fluctuaciones aleatorias.
A continuación, se presenta el análisis gráfico de las variables y la extracción de señales mediante la descomposición de las tres series de tiempo en estudio. Esto permitirá interpretar en detalle los patrones emergentes, perfeccionar la precisión de los pronósticos y generar insights estratégicos que fortalezcan la competitividad de Juan Valdez en el dinámico mercado del café.
Gráfico inicial de la variable Produccion de Cafe expresada en miles de sacos de 60 kg - Serie Original
La serie de tiempo objeto de análisis en este documento corresponde al nivel de producción de café, expresado en miles de sacos (60 kg). Entender su evolución es crucial, pues incide directamente en el rendimiento comercial y financiero de Juan Valdez. Las fluctuaciones en la oferta afectan la disponibilidad del grano y, por ende, los costos de materia prima, lo que exige ajustes en la estrategia de precios para preservar los márgenes de rentabilidad sin afectar el posicionamiento de la marca en los mercados nacionales e internacionales ni obstaculizar sus planes de expansión. Una disminución en la producción, por ejemplo, puede provocar escasez, elevar costos y restringir el crecimiento, mientras que una producción robusta favorece la estabilidad y el desarrollo. Además, la dependencia del café 100% colombiano (principal distintivo de la marca) expone a la empresa a riesgos climáticos y fitosanitarios que podrían deteriorar la calidad del grano en períodos de baja producción, obligándola a recurrir a materia prima de inferior calidad, ante una la limitada disponibilidad del insumo, y comprometiendo la percepción premium del producto.
Por ello, utilizar el análisis de series de tiempo y extraer señales (tendencia, estacionalidad y residuo) resulta crucial para detectar patrones, prever variaciones y optimizar la toma de decisiones estratégicas. Esto podría abarcar la formalización de contratos de compra de café a futuro, la adopción de coberturas financieras para mitigar la volatilidad de los costos o la definición de prioridades en la producción de líneas de productos específicas. Asimismo, ante posibles irrupciones en la oferta, es esencial robustecer los controles y desarrollar programas de aseguramiento de la calidad del producto. En el ámbito de la inversión y el crecimiento, sería necesario reconsiderar la priorización de la apertura de nuevas tiendas y la expansión hacia nuevos mercados, así como revisar y ajustar las estrategias de precios, marketing y ventas en los canales estratégicos que se deseen potenciar.
library(ggplot2)
library(plotly)
# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
Cafe$variable1 <- as.numeric(variable1_ts)
# Crear el gráfico
grafico_serie <- ggplot(Cafe, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(Cafe)), y = variable1)) + geom_line(color = "grey", linewidth = 0.4) + # Cambiado 'size' por 'linewidth'
geom_point(color = "black", size = 0.1) +
ggtitle("Produccion de Cafe: Serie original") +
xlab("Año") +
ylab("Miles de sacos (60kg)") +
theme_minimal()
ggplotly(grafico_serie)En la grafica se presenta la serie de tiempo mensual de la producción de café, expresada en miles de sacos (60 kg). Con el objetivo de comprender los factores que impulsan las variaciones observadas, se procede a descomponer la serie en sus componentes fundamentales: tendencia, estacionalidad y ruido. Esta metodología permite identificar de manera precisa posibles cambios estructurales en la serie.
Extracción señales de la variable Producción de café, expresada en miles de sacos de 60 kg.
# 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 Produccion de Cafe",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)La variable produccion de cafe, exhibe en su componente estacional, mayores volúmenes de producción hacia los últimos meses del año y niveles más bajos alrededor del quinto mes. Asimismo, el residuo se caracteriza por una alta volatilidad, evidenciándose al menos tres o cuatro choques atípicos, lo que podría complicar la selección del modelo de predicción óptimo, dado que los modelos ARIMA y SARIMA asumen un error de tipo ruido blanco. Por ello, es indispensable realizar un análisis detallado del error y aplicar las transformaciones pertinentes a la serie, con el fin de aprovechar al máximo las capacidades predictivas del modelo. Finalmente, se logra identificar con claridad el componente de tendencia, el cual resulta crucial para comprender el comportamiento estructural de la serie. En este componente, se observa inicialmente un período de tendencia creciente desde 2012 hasta finales de 2015, seguido de una fase de movimiento lateral. Posteriormente, se evidencia una tendencia decreciente durante los años 2022 y 2023, alcanzando un nivel mínimo durante el primer semestre del 2023, para finalmente registrar, aparentemente, un cambio hacia una tendencia creciente a finales de 2023.
Extracción señales Precio externo del cafe en centavos de dolar por libra
# 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 Precio del Cafe",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Al descomponer la serie del precio del café, se observa que los picos del componente estacional se alinean con los valles de la serie de producción, lo que sugiere una correlación entre los incrementos estacionales del precio y las disminuciones estacionales en la producción. No obstante, el análisis del componente de error revela periodos de considerable volatilidad, lo que destaca la necesidad de explorar otras variables que puedan influir en el comportamiento del precio antes de formular conclusiones definitivas. Finalmente, se resalta el cambio en la tendencia detectado a fines de 2023, evidenciado por un comportamiento ascendente progresivo que se mantiene hasta diciembre de 2024, fecha hasta la cual se disponen datos.
Extracción señales Exportacion de cafe en miles de dolares
# 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 exportaciones de cafe",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Al descomponer la serie de exportaciones de café y comparar sus patrones con los observados en las series anteriores, se evidencia un comportamiento estacional similar al de la producción de café, caracterizado por picos hacia finales de año y valles alrededor del quinto mes. Además, se identifican periodos de volatilidad pronunciada en el componente de error, consistente con lo observado en otras variables. Finalmente, la tendencia de las exportaciones se asemeja a la del precio del café, lo que sugiere que los incrementos en el precio podrían explicarse, al menos en parte, por un aumento en la demanda internacional de este producto.
Después de la descomposición temporal de cada variable, se extrae la variable ajustada por estacionalidad para graficarla junto con la serie original:
Se crea la variable Produccion de cafe ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]Se crea la variable Precio del cafe ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]Se crea la variable Exportacion de cafe ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]Ahora sí se pueden graficar las series originales versus las series ajustadas por estacionalidad
Gráfico serie original VS ajustada: Produccion de cafe
# 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("Produccion de Cafe: Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Produccion de cafe miles de sacos (60kg)") +
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)Al representar gráficamente la serie original junto con la serie ajustada por estacionalidad para la producción de café, se observa que, a pesar de que el ajuste atenúa el impacto de las fluctuaciones recurrentes, persiste un componente estacional notable. La presencia de un patrón estacional robusto tras el ajuste sugiere que los factores estacionales (como los ciclos climáticos y periodos de cosecha) tienen una influencia significativa en la dinámica de la serie. Esto podría indicar, además, la incidencia de otros factores exógenos, lo que justifica continuar la investigación para dilucidar con mayor precisión cómo el estudio de variables adicionales o la aplicación de modelos multivariables con variables exógenas puede enriquecer la explicación del comportamiento en la producción de café.
Gráfico serie original VS serie ajustada del precio del cafe
# Crear vector de fechas correctamente alineado con la serie
fechas_var2 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var2 <- ggplot() +
geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Variable 2:Serie Original vs Serie Ajustada del precio del cafe por Estacionalidad") +
xlab("AñoS") +
ylab("Centavos de dolar por libra") +
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)Al comparar la serie original con la serie ajustada por estacionalidad para el precio del café, se observa una notable similitud en sus comportamientos, lo que evidencia la ausencia de un patrón estacional definido. Mientras que la producción de café puede estar condicionada por ciclos climáticos y estacionales, en teoría, los precios se rigen por múltiples factores, como la oferta y la demanda global, los precios de futuros, la política económica y eventos macroeconómicos. En un escenario de fuerte estacionalidad, se esperaría identificar patrones cíclicos claros en la serie original; sin embargo, la mínima diferencia observada entre la línea que representa gráficamente la serie original del precio del café y la serie ajustada por estacionalidad indica que la variabilidad de estos precios se debe mayormente a factores aleatorios o estructurales y no a un componente estacional. Esto resalta la relevancia de investigar cómo otras variables exógenas pueden influir en el comportamiento del precio del café, ya que estos podrían estar determinados por elementos de mayor volatilidad y alcance global.
Gráfico serie original VS serie ajustada exportacion de cafe
# 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 exportacion de cafe") +
ggtitle("Variable 3:Serie Original vs Serie Ajustada exportacion de cafe por Estacionalidad") +
xlab("Años") +
ylab("Miles de dolares") +
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)Similar a lo observado en la producción de café, la estacionalidad es un componente robusto en la serie de exportaciones de café, lo que indica que los patrones estacionales están profundamente arraigados en el comportamiento de ambas variables. La persistencia de un componente estacional significativo tras el ajuste sugiere que factores como los ciclos climáticos, los periodos de cosecha y los patrones de demanda ejercen una influencia considerable. Este hallazgo podría justificar la incorporación de variables exógenas y la aplicación de modelos multivariantes para capturar con mayor precisión la dinámica subyacente en la producción y exportación de café. Asimismo, se recomienda considerar el uso de modelos especializados en el manejo de la estacionalidad y en la incorporación de factores externos, lo que contribuirá a mejorar la estimación de las variables en estudio.
Ahora graficamos las series originales vs tendencia
El análisis gráfico comparativo entre las series originales de producción, precio y exportaciones de café y sus respectivas tendencias permite identificar de forma precisa los cambios estructurales subyacentes en cada variable. La extracción de la tendencia resulta fundamental para prever escenarios futuros y anticipar crisis u oportunidades emergentes en el sector cafetalero. Estos hallazgos ofrecen una base sólida para el desarrollo de estrategias financieras y comerciales, lo que puede orientar a la empresa Juan Valdez en la formulación de políticas adaptadas a las dinámicas del mercado y a la evolución de la industria.
Primero se debe obtener la tendencia de cada variable y luego graficarla
Tendencia Produccion de Cafe
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("Produccion de Cafe: Serie Original vs Tendencia") +
xlab("Años") +
ylab("Miles de sacos (60kg)") +
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)El análisis revela un período de tendencia creciente sostenido desde 2012 hasta finales de 2015, seguido de una fase de comportamiento lateral. Posteriormente, se identifica una tendencia a la baja durante 2022 y 2023, alcanzando valores mínimos en el primer semestre de 2023, para finalmente registrar, aparentemente, un giro hacia una tendencia ascendente a fines de 2023.
Con base en estos hallazgos, se recomienda implementar estrategias financieras que incluyan la diversificación de las fuentes de financiación, lo que permitiría capitalizar la tendencia creciente en la producción y facilitar el impulso de proyectos de expansión en nuevos mercados, el lanzamiento de nuevos productos o el desarrollo de canales estratégicos. Paralelamente, es fundamental adoptar instrumentos de cobertura que mitiguen los riesgos asociados a la volatilidad del mercado, asegurando así una estructura de costos estable y protegiendo el margen de rentabilidad.
En el ámbito comercial, resulta crucial optimizar la cadena de suministro y ajustar la política de precios para aprovechar las señales de recuperación del mercado. Por último, en términos de marketing, se sugiere reforzar el posicionamiento de la marca Juan Valdez mediante campañas que destaquen la calidad y el origen premium del café, orientadas tanto a captar nuevos segmentos de mercado como a fidelizar a la clientela existente.
Tendencia Precio del Cafe
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("Precio del Cafe: Serie Original vs Tendencia") +
xlab("Años") +
ylab("Centavos de dolar por libra") +
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)El análisis del gráfico que compara la tendencia del precio con la serie original revela un giro hacia una tendencia ascendente a fines de 2023, que se mantiene de forma progresiva hasta diciembre de 2024, posiblemente impulsado por un aumento en la demanda internacional. Es esencial confirmar si este cambio constituye un verdadero cambio estructural en la serie, lo que permitiría diseñar una estrategia de precios que capture el valor percibido por el mercado. Simultáneamente, se recomienda implementar un sistema de precios dinámico, adaptable a los ciclos estacionales y a las variaciones en la oferta, para aprovechar los picos de demanda y mitigar el impacto de la escasez de materia prima. Además, la adopción de estrategias de cobertura mediante el uso de derivados financieros ayudará a proteger los márgenes de rentabilidad frente a fluctuaciones inesperadas. En caso de que se verifique que el aumento de precios se debe a un cambio estructural, sería conveniente fortalecer el posicionamiento premium de la marca mediante campañas que resalten la calidad y el origen del café, orientadas a captar nuevos segmentos de mercado, especialmente a nivel internacional.
Tendencia Exportacion de cafe
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("Exportacion de cafe: Serie Original vs Tendencia") +
xlab("Año") +
ylab("Miles de dolares") +
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)Dado que la dinámica de las tendencias del precio del café y las exportaciones presenta un comportamiento notablemente similar a partir del 2020, las estrategias comerciales y financieras propuestas resultan compatibles para las dos variables analizadas.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia:
Tasa de crecimiento de la serie de tendencia y original para la Produccion del cafe
#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] 144print(length(tasa_crecimiento_var1))## [1] 144print(length(tasa_tendencia_var1))## [1] 144*Gráfico variable original y tendencia Produccion de cafe: 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("Produccion de Cafe: Tasa de crecimiento anual % de la serie Original y la tendencia") +
xlab("Años") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)El aumento sostenido en la tasa de crecimiento de la producción de café desde finales del 2023, indica que esta variable se encuentra en un periodo de expansión. Esto representa una gran oportunidad para Juan Valdez, y esta muy alineado con su plan estrategico de expansion en mercados internacionales.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: Precio del Cafe
#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] 144print(length(tasa_crecimiento_var2))## [1] 144print(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("Precio del Cafe: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
xlab("Años") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)La sostenida tendencia ascendente en la tasa de crecimiento del precio del café, observada desde mediados de 2024, respalda de forma contundente las iniciativas de crecimiento y expansión de la empresa Juan Valdez, al reflejar una demanda robusta que inspira confianza en el mercado.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: Exportacion de cafe
#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] 144print(length(tasa_crecimiento_var3))## [1] 144print(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("Exportacion de cafe: Tasa de crecimiento anual % de la serie Original y la tendencia") +
xlab("Años") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)El análisis integral de las series de producción, precio y exportaciones de café demuestra que Juan Valdez opera en un entorno dinámico, caracterizado por una expansión en la producción y una tendencia ascendente en los precios desde mediados de 2024, acompañada por una creciente demanda internacional durante el mismo periodo. No obstante, la persistente volatilidad en los residuos evidencia la influencia de factores exógenos, lo que subraya la necesidad de implementar modelos predictivos multivariantes para detectar otras señales del mercado. En este contexto, resulta esencial que Juan Valdez establezca estrategias de cobertura financiera que aseguren una estructura de costos estable. Se recomienda adoptar una política de precios dinámica, diseñada para proteger los márgenes de rentabilidad y adaptable tanto a los ciclos estacionales como a las variaciones en la oferta, además de optimizar la cadena de suministro para mejorar la eficiencia operativa y garantizar la calidad del producto. Asimismo, diversificar las fuentes de financiamiento se revela como una estrategia clave para respaldar proyectos de expansión, en especial en mercados internacionales, y fomentar el desarrollo de nuevas líneas de productos, aprovechando el actual momentum del mercado. Paralelamente, es crucial reforzar el posicionamiento premium de la marca a nivel global mediante campañas de marketing que resalten la calidad y el origen 100% colombiano del café, con el fin de captar nuevos segmentos de mercado y fidelizar a la clientela existente. Estas iniciativas permitirán a Juan Valdez anticipar cambios en el entorno, gestionar la volatilidad y capitalizar oportunidades de crecimiento sostenible, consolidándose como una marca de referencia en un entorno global competitivo.
Modelo ARIMA para pronosticar la producción de café
A continuación, se presenta la metodología adoptada para seleccionar un modelo ARIMA idóneo para pronosticar la producción de café, siguiendo el enfoque Box-Jenkins. Inicialmente, se realizó un análisis del comportamiento histórico de la serie temporal, identificando patrones, tendencias y componentes estacionales. Ahora, en esta segunda etapa, se aplicará el test ADF para evaluar la estacionariedad de la variable “producción de café” y determinar las transformaciones necesarias para estabilizar la serie original. Una vez definida la serie, ya sea en su forma transformada o original, se procederá a evaluar las funciones de autocorrelación y autocorrelación parcial para establecer los parámetros iniciales del modelo. La etapa final abarcará la estimación y validación del modelo, asegurando que los residuos se comporten como ruido blanco. Finalmente, se seleccionará el modelo que ofrezca la mayor precisión y confiabilidad en el pronóstico, sirviendo como base para respaldar decisiones estratégicas en la gestión de la producción de café.
Desglose del nombre ARIMA
- Autoregressive o autorregresivo (AR) → Usa valores pasados para predecir el futuro
- Integrated (I) → Ajusta tendencias en los datos (sin patrones cambiantes en el tiempo). 3.Moving Average o media móvil (MA)- → Suaviza fluctuaciones aleatorias a partir de errores pasados.
ARIMA combina estos tres elementos para crear una predicción más precisa.
Metodología Box-Jenkins para aplicar un modelo ARIMA
La metodología Box-Jenkins es un enfoque sistemático para construir modelos ARIMA con el objetivo de analizar y pronosticar series de tiempo. Fue desarrollada por George Box y Gwilym Jenkins y se basa en cuatro etapas clave:
1️⃣ Identificación 2️⃣ Estimación 3️⃣ Validación 4️⃣ Pronóstico
Se usa especialmente cuando se quiere encontrar el modelo ARIMA más adecuado para una serie de tiempo.
Antes de empezar a aplicar la metodología BOX-JENKINS, lo ideal es dividir el conjunto de datos de prueba y entrenamiento
✅ Para entrenar el modelo con datos históricos sin usar información futura. ✅ Para evaluar la precisión del modelo comparando sus predicciones con los datos reales de prueba.
💡 Esta división es clave en modelos predictivos para evitar sobreajuste y evaluar el rendimiento en datos no vistos.
División en conjunto de entrenamiento y prueba para la variable de la producción de café que es la elegida para pronosticar
El código siguiente divide una serie temporal 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.
NOTA: Debido a la elevada volatilidad observada en la variable producción de café, se optó por trabajar con su logaritmo. Esta transformación estabiliza la varianza y convierte relaciones multiplicativas en aditivas, lo que facilita el análisis. Además, el uso del logaritmo mitiga el impacto de valores extremos, contribuyendo en última instancia a la formulación de un modelo más robusto y preciso.
log_variable1_ts <- log(variable1_ts)# Esta división idealmente podria ser 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(log_variable1_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(log_variable1_ts, end = c(2024, 9)) # Entrenamiento hasta septiembre 2024
test_ts <- window(log_variable1_ts, start = c(2024, 10)) # Prueba inicia desde oct2024Paso 1: Identificación del modelo
Identificar estacionariedad
En el análisis de series de tiempo, una serie es estacionaria si su comportamiento es constante a lo largo del tiempo, es decir:
✅ Su promedio no cambia con el tiempo. ✅ Su variabilidad (qué tanto fluctúa) se mantiene estable. ✅ Su relación con valores pasados es siempre la misma.
¿Cómo saber si una serie es estacionaria?
1️⃣ Observando un gráfico Si el gráfico de la serie muestra una tendencia creciente o decreciente, o si las variaciones se hacen más grandes con el tiempo, la serie probablemente no es estacionaria.
2️⃣ Usando la prueba de Dickey-Fuller Aumentada (ADF) Es una prueba estadística que nos dice si la serie tiene una tendencia fuerte. Si el p-valor de la prueba es menor a 0.05, significa que la serie es estacionaria.
¿Cómo hacer una serie estacionaria?
Si encontramos que la serie no es estacionaria, podemos transformarla para que lo sea:
✅ Diferenciación: Restamos cada valor con su valor anterior. Esto elimina tendencias crecientes o decrecientes. ✅ Tomar el logaritmo: Si la variabilidad crece con el tiempo, aplicar un logaritmo estabiliza la varianza. ✅ Eliminar tendencias o ajustar estacionalidad: Si hay patrones repetitivos, podemos restarlos o modelarlos por separado.
Test de Dickey-Fuller
El test de Dickey-Fuller aumentado (ADF) se usa para verificar si una serie temporal es estacionaria, es decir, si sus propiedades estadísticas (media y varianza) permanecen constantes en el tiempo.
HO: Serie no estacionaria HI: Serie estacionaria
¿Qué significa el p-valor?
Si el p-valor es bajo (< 0.05) → Rechazamos la hipótesis nula y concluimos que la serie es estacionaria. Si el p-valor es alto (> 0.05) → No podemos rechazar la hipótesis nula, lo que indica que la serie no es estacionaria.
A continuación se aplica el test ADF para validar estacionariedad en el conjunto de entrenamiento de la variable logaritmo de la producción de café, la cual es la elegida para pronosticar la producción de café
library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable 1 (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test##
## Augmented Dickey-Fuller Test
##
## data: train_ts
## Dickey-Fuller = -3.0901, Lag order = 5, p-value = 0.1221
## alternative hypothesis: stationaryEl test ADF en la variable logaritmo de la producción de café, arrojó un p-value igual a 0.1221, este valor es mayor a 0.05, por tanto la serie es no estacionaria. De ese modo se debe ejecutar el código siguiente para diferenciar una vez la variable logaritmo de la producción de café y luego volver a aplicar el test ADF a esa serie diferenciada una vez:
#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable 1 , una sola vez:
train_diff <- diff(train_ts, differences = 1) Diferenciación en niveles de la variable logaritmo de la producción de café
A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable logaritmo de la producción de café para ver graficamente el cambio o ajuste:
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "blue") +
ggtitle("Produccion de cafe:Serie Original") +
xlab("Año") + ylab("Produccion de cafe")
ggplotly(p2) # Convertir en gráfico interactivo # Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], variable1_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable1_Diff)) +
geom_line(color = "red") +
ggtitle("Log Produccion de cafe: Serie Estacionaria (Una diferenciación)") +
xlab("Año") + ylab("Log Produccion de cafe 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoAl comparar ambas representaciones, se evidencia la eliminación de tendencias y, dado que la variable está expresada en logaritmos, se observa una notable reducción en la volatilidad de la serie.
Ahora probamos estacionariedad en la serie diferenciada (en logaritmo)
# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
adf_train_diff <- adf.test(train_diff)
print(adf_train_diff)##
## Augmented Dickey-Fuller Test
##
## data: train_diff
## Dickey-Fuller = -7.3287, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEl test ADF en la variable diferenciada una vez, de la variable logaritmo de la producción de café, arrojó un p-value igual a 0.01, este valor es manor a 0.05, por tanto la serie ES ESTACIONARIA. Se procede a trabajar entonces en un modelo log-ARIMA con una diferenciación (d=1).
Identificación manual de p y q
ACF (Autocorrelation Function)
Muestra la correlación de la serie con sus rezagos.
Ayuda a determinar el parámetro q en un modelo ARIMA(p, d, q).
PACF (Partial Autocorrelation Function)
Muestra la correlación parcial entre la serie y un rezago específico, eliminando el efecto de rezagos intermedios.
Ayuda a determinar el parámetro p en un modelo ARIMA(p, d, q).
Recordemos El eje X representa los rezagos (lags). El eje Y muestra la autocorrelación parcial en cada rezago. Las líneas azules punteadas son los intervalos de confianza (aproximadamente 95%). Si una barra sobrepasa estos límites, indica una autocorrelación significativa. Si las barras caen dentro de los límites, no son significativamente diferentes de cero
En el código siguiente se crean los correlogramas para determinar los posibles valores que puedeo tomar el parámetro p y q:
library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")
ggplotly(acf_plot)ggplotly(pacf_plot)Interpretación correlogramas
Se puede observar que los valores que podrian tomar p y q serian:
p=2 p=4* (P optimo=4) q=2 q=4* q=6 (q óptimo=4)
El modelo óptimo para esta variable seria (4,1,4)
AIC y BIC: Se usan para comparar modelos; cuanto más bajos, mejor. Las métricas AIC y BIC son fundamentales para la comparación de modelos, ya que penalizan la complejidad y favorecen aquellos que logran un buen ajuste con menos parámetros; es decir, mientras más bajos sean estos valores, mejor es la relación entre ajuste y parsimonia del modelo. Durante el proceso de iteración se evaluó si existía un modelo superior al ARIMA (4,1,4). Tras comparar los valores de AIC, AICc y BIC, se constató que el modelo ARIMA (0,1,4) arrojaba cifras inferiores en todas las métricas, lo que indica un desempeño más óptimo. En consecuencia, se seleccionó el modelo ARIMA (0,1,4) como el modelo más óptimo según estos criterios de selección.
Nuevo modelo optimo segun las metricas de seleccion de modelos AIC, AICc y BIC—>ARIMA (0,1,4)
Paso 2. Estimación manual del modelo
Métricas de evaluación
Mean Absolute Error (MAE)
Representa el error absoluto promedio entre las predicciones del modelo y los valores reales.
Root Mean Squared Error (RMSE) = 0.1588512
El RMSE de 0.1588512 señala que el error cuadrático medio, medido en las unidades de la variable (en la escala logarítmica), es muy bajo, lo que confirma que la magnitud absoluta del error es reducida.
Mean Absolute Percentage Error (MAPE) = 1.86%
Expresa el error en términos relativos, como porcentaje del valor real.
Interpretación: Un MAPE de 1.857125 indica que, en promedio, las predicciones del modelo se desvían aproximadamente un 1.86% respecto a los valores reales. Esto sugiere que el modelo presenta una alta precisión en sus pronósticos, ya que el error porcentual absoluto medio es muy bajo.
Regla general: MAPE < 10% → Muy buen modelo ✅ 10%-20% → Modelo aceptable 👍 20%-50% → Modelo pobre ⚠️ 50% → Modelo muy malo ❌
En este caso, un MAPE de 1.86% sugiere un muy buen modelo para pronostico.
Estimación del modelo identificado (0,1,4)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(0,1,4)) #Se va a estimar un modelo Arima de orden (0,1,4)
summary(manual_arima_model)## Series: train_ts
## ARIMA(0,1,4)
##
## Coefficients:
## ma1 ma2 ma3 ma4
## -0.2238 -0.4525 0.0389 -0.1309
## s.e. 0.0848 0.1085 0.0772 0.1070
##
## sigma^2 = 0.02609: log likelihood = 63.01
## AIC=-116.02 AICc=-115.61 BIC=-100.9
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01630767 0.1588512 0.1285159 0.2002336 1.857125 0.7809862
## ACF1
## Training set -0.02960454Significancia de coefientes Dado que este en este modelo no se incluyen términos autorregresivos (p = 0), no se incorporan coeficientes para estimar la insidencia del componente autorregresivo. Por el contrario, se incorporan cuatro términos de media móvil (q = 4) para modelar la estructura del error. Bajo una significancia del 11%, todos los coeficientes son significativos.
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.223750 0.084773 -2.6394 0.008305 **
## ma2 -0.452483 0.108455 -4.1721 3.018e-05 ***
## ma3 0.038877 0.077245 0.5033 0.614757
## ma4 -0.130949 0.107016 -1.2236 0.221091
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación significancia coeficientes
Se puede observar que MA1 y MA2, son altamente significativos, lo que significa que estos coeficientes de MA tienen un impacto importante en el modelo, en contraste, MA3 y MA4 no son significativos. Como al menos dos de los cuatro componentes son significativo, entonces continuo con la validación de residuos del modelo.
Paso 3. Validación de residuos del modelo estimado manual
La validación de residuos es crucial para determinar si el modelo ARIMA es adecuado o si necesita mejoras. El objetivo es verificar que los residuos (errores de predicción) se comporten como ruido blanco, es decir, sin patrones detectables.
1. Serie de residuos (gráfico superior)
Muestra cómo se comportan los errores a lo largo del tiempo. Idealmente, deberían oscilar alrededor de cero sin tendencias evidentes ni grandes acumulaciones de error.
Posible problema identificado: Se aprecian caídas pronunciadas en torno al 2020 y al 2021, lo que podría señalar la ocurrencia de cambios estructurales en los datos, posiblemente atribuibles a eventos como la pandemia y el paro nacional. Adicionalmente, el gráfico del error evidencia una alta volatilidad, lo cual indica la presencia de una variabilidad significativa que el modelo no logra explicar, sugiriendo la influencia de factores estructurales o de eventos atípicos. Esta volatilidad elevada puede comprometer la precisión de las predicciones, lo que resalta la necesidad de revisar y, en su caso, ajustar el modelo para captar de manera más adecuada la dinámica subyacente de la serie. Alternativamente, se podría considerar la adopción de un modelo multivariante que integre variables exógenas, permitiendo así explicar mejor el comportamiento de la producción de café, aspecto que el modelo univariante actual no logra abordar completamente.
Función de Autocorrelación (gráfico inferior izquierdo, ACF de residuos)
Si el modelo es adecuado, los residuos no deben mostrar correlaciones significativas en el tiempo.
Interpretación: Al analizar el gráfico de la función de autocorrelación (ACF) se observa que en varias ocasiones los componentes de error se ubican fuera de los intervalos de confianza, esto indica que los errores presentan autocorrelación significativa. En otras palabras, los errores no son independientes, lo que sugiere que el modelo no ha capturado completamente la estructura temporal de la serie y que existen patrones sistemáticos en los residuos que podrían ser aprovechados para mejorar la especificación del modelo.
Histograma de residuos con ajuste normal (gráfico inferior derecho)
Sirve para verificar si los errores siguen una distribución normal, lo cual es un supuesto clave en ARIMA. Interpretación: La curva roja representa la distribución normal teórica. El análisis del histograma de la distribución de errores revela que, si bien los residuos se aproximan a una distribución normal, se evidencia un leve sesgo hacia la izquierda y la presencia de valores extremos, manifestados en colas más densas de lo esperado. Esto sugiere la existencia de eventos atípicos o de datos que no son completamente capturados por el modelo actual, lo que podría afectar la precisión de las predicciones. Dichos hallazgos indican la necesidad de explorar métodos robustos o técnicas de tratamiento de outliers, para mejorar el ajuste y la capacidad explicativa del modelo.
Conclusión y acciones recomendadas
✅ A pesar de que el modelo ARIMA(0,1,4) presentó resultados prometedores en términos de RMSE y MAPE, no se considera óptimo, ya que el análisis de los residuos revela que estos no se comportan como ruido blanco. Esto sugiere la presencia de eventos atípicos o información que el modelo no logra capturar adecuadamente. Por ello, se recomienda explorar métodos para el tratamiento de los errores o considerar la aplicación de modelos alternativos, especialmente aquellos diseñados para series con fuertes componentes estacionales. Asimismo, es pertinente investigar la viabilidad de modelos multivariantes que puedan explicar con mayor precisión la producción de café.
Recordar: El supuesto de normalidad significa que los errores o residuos de un modelo deben seguir una distribución normal (o “campana de Gauss”). Si los errores son normales, podemos hacer predicciones más confiables y usar ciertas pruebas estadísticas que asumen esta propiedad.
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,4)
## Q* = 76.813, df = 20, p-value = 1.355e-08
##
## Model df: 4. Total lags used: 24Paso 4. Pronóstico (modelo manual)
Al comparar el pronostico del modelo manual vs. los valores observados, se evidencia que a pesar de que se captura el cambio en la tendencia, el pronostico no sigue los movimientos de la tendecia y adicionalmente presenta un alto nivel de subestimacion.
Como se planteaba al principio del documento, esta serie al presentar patrones estacionales fuertes y no estar bien capturados, el modelo puede fallar en prever las fluctuaciones. En ese caso se sugiere, nuevamente, trabajar con otro modelo como SARIMA, o modelo multivariante que incorpore variables exogenas que permitan pronosticar de manera mas precisa la produccion del cafe.
Pronóstico en el test de prueba (oct, nov y dic 2024) y gráfico
#Aquí se crea el pronóstico con el modelo ARIMA manual y se guarda en una nueva riable u objeto "manual_forecast"
manual_forecast <- forecast(manual_arima_model, h = length(test_ts)) #Se generan tantos pronósticos como valores tenga el conjunto de prueba (test_ts).
# Crear dataframe para gráfico interactivo del pronóstico manual
manual_forecast_data <- data.frame(Tiempo = time(manual_forecast$mean), ## Extrae las fechas del pronóstico
Pronostico = as.numeric(manual_forecast$mean), ## Valores pronosticados
Observado = as.numeric(test_ts)) ## Valores reales de la serie
# Graficar pronóstico manual junto con los valores observados reales
p_manual <- ggplot(manual_forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico Manual")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Log.Produccion Cafe:Pronóstico Manual vs Observado") +
xlab("Año") + ylab("Produccion Cafe miles de sacos (kg)")
ggplotly(p_manual)Métricas de evaluación del modelo manual dentro del periodo de prueba (oct,nov y dic2024
# Calcular métricas de evaluación del modelo manual
mae_manual <- mean(abs(manual_forecast$mean - test_ts), na.rm = TRUE)
rmse_manual <- sqrt(mean((manual_forecast$mean - test_ts)^2, na.rm = TRUE))
# Mostrar métricas de evaluación del modelo manual
cat("MAE Manual: ", mae_manual, "\n")## MAE Manual: 0.4370279cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 0.4652763A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic2024
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_manual <- forecast(manual_arima_model, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table_manual <- data.frame(
Tiempo = time(arima_forecast_manual$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast_manual$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table_manual)## Tiempo Observado Pronosticado
## 1 2024.750 7.199794 6.988514
## 2 2024.833 7.473870 6.927354
## 3 2024.917 7.494558 6.941270En la tabla se evidencia la subestimacion en cada uno de los periodos.
Ahora pronosticamos fuera del periodo de análisis: Enero 2025
Es decir, le sumamos al periodo de prueba (oct,nov,dic2024) una observación más (enero 2025). 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_manual <- forecast(manual_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
next_month_forecast_manual <- data.frame(
Tiempo = time(next_forecast_manual$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast_manual$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast_manual)## Tiempo Pronostico
## 1 2024.750 6.988514
## 2 2024.833 6.927354
## 3 2024.917 6.941270
## 4 2025.000 6.928997# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_manual, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2025: 2025 = 6.92899686122672"Se confirma que al menos, el modelo captura el cambio estacional que sucede en el primer mes del año.
Otra forma para calcular un valor futuro (fuera de muestra)-Modelo manual, es decir, en caso de que no se haga la dviisón inicial de conjunto de entrenamiento y prueba
Si no hay conjunto de prueba o test y solo quieres el siguiente punto u observación, el código a usar seria algo asi:
# Pronosticar octubre de 2024
future_forecast_manual <- forecast(manual_arima_model, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024 <- future_forecast_manual$mean[1]
print(paste("Pronóstico para oct 2024:", forecast_oct2024))## [1] "Pronóstico para oct 2024: 6.98851378652133"Modelo ARIMA automático
Usa la función auto.arima() de forecast en R para seleccionar automáticamente los mejores parámetros (p,d,q).
✅ Ventajas:
✔ Optimización automática: Encuentra los valores óptimos de ARIMA sin intervención manual. ✔ Ahorra tiempo: Útil cuando hay muchas series a modelar. ✔ Evita sesgo humano: Reduce el riesgo de elegir un modelo incorrecto por falta de experiencia. ✔ Incluye corrección por estacionalidad si se usa con seasonal = TRUE. ✔ Suele funcionar bien en la mayoría de los casos, ya que usa criterios como AIC/BIC para optimizar.
❌ Desventajas: ❌ Puede no ser el mejor modelo posible, ya que depende del criterio de selección. ❌ Menos interpretabilidad: No siempre es claro por qué eligió ciertos parámetros. ❌ Puede ignorar conocimiento experto sobre la serie o factores externos.
Identificación automática del modelo ARIMA
El modelo automático identificado es (0,1,2). Si se compara el AIC o BIC de este modelo frente el modelo manual (0,1,4), se obtiene un valor más bajo en esta métricas obtenidas al estimar el modelo manualmente. Por tanto, el modelo ARIMA (0,1,4) sigue siendo preferieble.
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(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.2315 -0.5336
## s.e. 0.0761 0.0747
##
## sigma^2 = 0.0261: log likelihood = 61.92
## AIC=-117.84 AICc=-117.68 BIC=-108.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01652494 0.159976 0.1294284 0.2027561 1.870519 0.7865316
## ACF1
## Training set -0.0459933Significancia de coeficientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.23152 0.07607 -3.0435 0.002338 **
## ma2 -0.53356 0.07468 -7.1446 9.025e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Ajuste del modelo ARIMA(4,1,2) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts,
order = c(0, 1, 2)) # Especificamos directamente (p=4, d=1, q=2)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.2315 -0.5336
## s.e. 0.0761 0.0747
##
## sigma^2 = 0.0261: log likelihood = 61.92
## AIC=-117.84 AICc=-117.68 BIC=-108.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.01652494 0.159976 0.1294284 0.2027561 1.870519 0.7865316
## ACF1
## Training set -0.0459933# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_auto) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)
## Q* = 77.692, df = 22, p-value = 3.869e-08
##
## Model df: 2. Total lags used: 24Pronóstico modelo ARIMA automático (0,1,2)
# Generar pronóstico para el conjunto de prueba
forecast_arima_auto <- forecast(darima_auto, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data_auto <- data.frame(Tiempo = time(forecast_arima_auto$mean),
Pronostico = as.numeric(forecast_arima_auto$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4auto <- ggplot(forecast_data_auto, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Pronóstico vs Observado") +
xlab("Tiempo") + ylab("variable1")
ggplotly(p4auto) # Convertir el gráfico en interactivoEl modelo automático (0,1,2) parece pronosticar peor dentro de prueba, ya que al menos el modelo (0,1,4) logra generar un ligero aumento a partir de los ultimos meses de año, mientras que el automatico muestra mas bien un movimiento lateral. Ambos modelos presentan un problema de subestimacion.
# 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 7.199794 6.983953
## 2 2024.833 7.473870 6.936842
## 3 2024.917 7.494558 6.936842Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueb auna observación más. Es decir, se estan pronosticando 4 observaciones o trimestres.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo trimestre
next_month_forecast_auto <- data.frame(
Tiempo = time(next_forecast_auto$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast_auto$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast_auto)## Tiempo Pronostico
## 1 2024.750 6.983953
## 2 2024.833 6.936842
## 3 2024.917 6.936842
## 4 2025.000 6.936842# 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 = 6.93684232877514"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar el mes octubre 2024
future_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024_auto <- future_forecast_auto$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_auto))## [1] "Pronóstico para octubre 2024: 6.98395279279621"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, como la produccion de cafe.
Un SARIMA (Seasonal ARIMA) es un modelo que combina:
✅ Un modelo ARIMA tradicional para capturar relaciones entre valores pasados y errores. ✅ Un componente estacional, útil cuando los datos muestran patrones repetitivos en el tiempo (por ejemplo, ventas trimestrales, datos mensuales de temperatura, etc.).
Básicamente, es como un ARIMA mejorado que puede manejar ciclos estacionales.
En este caso, el modelo detectó un patrón repetitivo cada 12 períodos, por lo que usa componentes autoregresivos y de media móvil para modelarlo
El modelo ajustado es un SARIMA(1,1,2)(2,0,0)[12], lo que significa:
(1,1,2): Parte ARIMA no estacional: 1 términos autorregresivos (AR). 1 diferenciación (d), lo que indica que la serie fue diferenciada una vez para hacerla estacionaria. 2 término de media móvil (MA).
(2,0,0)[12]: Parte estacional con periodicidad 12 (mensual): 2 término autorregresivo estacional (SAR). 0 diferenciaciones estacionales. 0 términos de media móvil estacionales (SMA).
# Identificación automática del mejor modelo ARIMA
auto_arima_model <- auto.arima(train_ts) # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)## Series: train_ts
## ARIMA(1,1,2)(2,0,0)[12]
##
## Coefficients:
## ar1 ma1 ma2 sar1 sar2
## -0.3111 -0.1201 -0.5631 0.3021 0.3371
## s.e. 0.1648 0.1338 0.0928 0.0786 0.0826
##
## sigma^2 = 0.0192: log likelihood = 84.03
## AIC=-156.06 AICc=-155.48 BIC=-137.91A continuación, se crea el objeto darima para lueg poder graficar los valores reales y observados:
# Cargar el paquete necesario
library(forecast)
# Ajustar el modelo SARIMA(1,1,2)(2,0,0)[12]
darima <- Arima(train_ts,
order = c(1, 1, 2), # (p,d,q) -> (0,1,1)
seasonal = list(order = c(2, 0, 0), # (P,D,Q) -> (1,0,0)
period = 12)) # Periodicidad estacional de 12 meses
# Mostrar resumen del modelo ajustado
summary(darima)## Series: train_ts
## ARIMA(1,1,2)(2,0,0)[12]
##
## Coefficients:
## ar1 ma1 ma2 sar1 sar2
## -0.3111 -0.1201 -0.5631 0.3021 0.3371
## s.e. 0.1648 0.1338 0.0928 0.0786 0.0826
##
## sigma^2 = 0.0192: log likelihood = 84.03
## AIC=-156.06 AICc=-155.48 BIC=-137.91
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.006896087 0.1358067 0.1069678 0.07348956 1.548463 0.6500396
## ACF1
## Training set 0.0006052778En el correlograma de residuos siguiente se observa que, mejora la correlación de los residuos frente a lso dos modelos anteriores (ver correlograma ACF donde los errores pasados se encuentran por dentro del intervalo de confianza).Sin embargo, no se logra reducir la volatilidad en el error debido a la presencia de datos atipicos. Finalmente, al comparar los valores reales vs. los pronosticados se determina una poca coincidencia, por lo que sigue funcionando mejor el modelo automatico (0,1,4)
Adicionalmente, al comparar las metricas AIC y BIC, el ARIMA (0,1,4) sigue siendo preferible.
# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2)(2,0,0)[12]
## Q* = 20.031, df = 19, p-value = 0.3927
##
## Model df: 5. Total lags used: 24Pronóstico con el modelo SARIMA
# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean),
Pronostico = as.numeric(forecast_arima$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Pronóstico vs Observado") +
xlab("Tiempo") + ylab("Unidad Variable 1")
ggplotly(p4) # Convertir el gráfico en interactivo# 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 7.199794 7.060928
## 2 2024.833 7.473870 7.097913
## 3 2024.917 7.494558 7.073478Ahora pronosticamos fuera del periodo de análisis
El problema de subestimacion persiste.
# 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 7.060928
## 2 2024.833 7.097913
## 3 2024.917 7.073478
## 4 2025.000 6.954323# 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 = 6.954323351684"Otra forma para calcular un valor futuro (fuera de muestra)
# Identificación automática del mejor modelo ARIMA
auto_arima_model <- auto.arima(train_ts) # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)## Series: train_ts
## ARIMA(1,1,2)(2,0,0)[12]
##
## Coefficients:
## ar1 ma1 ma2 sar1 sar2
## -0.3111 -0.1201 -0.5631 0.3021 0.3371
## s.e. 0.1648 0.1338 0.0928 0.0786 0.0826
##
## sigma^2 = 0.0192: log likelihood = 84.03
## AIC=-156.06 AICc=-155.48 BIC=-137.91# Pronosticar octubre 2024
future_forecast_sarima <- forecast(auto_arima_model, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024_sarima <- future_forecast_sarima$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_sarima))## [1] "Pronóstico para octubre 2024: 7.06092841094476"Conclusión:
El análisis comparativo de diversos modelos de predicción para la producción de café ha permitido constatar que el modelo manual ARIMA(0,1,4), presenta un rendimiento superior de acuerdo con los criterios de información de Akaike (AIC) y Bayesiano (BIC). No obstante, un estudio minucioso de los residuos pone de manifiesto una autocorrelación significativa y persistente, lo cual evidencia que estos no se comportan como un proceso de ruido blanco. Esta situación indica que el modelo no captura de forma adecuada la estructura inherente a la serie temporal, afectando así la validez y confiabilidad de los pronósticos obtenidos.
Aunque el modelo ARIMA(0,1,4) optimiza ciertos índices de bondad de ajuste frente a los otros modelos analizados, su desempeño se ve restringido por deficiencias en la modelización de los errores y la presencia de patrones estacionales fuertes y de datos atípicos. Por lo tanto, se vuelve imprescindible el desarrollo de modelos más robustos que aborden estas limitaciones.
En este contexto, se recomienda la exploración de enfoques multivariantes que integren variables exógenas pertinentes, lo cual podría facilitar una estimación más precisa y confiable de la producción de café. Además, la incorporación de técnicas avanzadas para la identificación y tratamiento de valores atípicos, junto con estrategias más sofisticadas para el manejo de la estacionalidad, podría incrementar sustancialmente la precisión de los pronósticos. Estos avances contribuirían no solo a perfeccionar el proceso de modelización, sino también a establecer un marco analítico más sólido para la toma de decisiones estratégicas basadas en un análisis estadístico riguroso.
By Pablo Lopez y Natalya Lopez.