EXTRACCIÓN DE SEÑALES Y MODELO ARIMA: EXPORTACIONES DE CAFÉ
Susan Jimenez y Patrick Simon
2026-04-24
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)
data_col <- read_excel("C:/Users/HONOR/Downloads/9NO_INTER_CALI/CASO2/Base Caso2.xlsx",
col_types = c("date", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric"))
View(data_col)PASO INDISPENSABLE: Declarar la (s) variable (s) como serie (s) temporal (es):
Variable 1
# Convertir/declarar variable 1 en serie de tiempo mensual
variable1_ts <- ts(data_col$XCAF, start = c(2012, 1), frequency = 12)
variable1_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
## 2025 474078.6 434404.7 550164.4 483948.2 427844.6 419873.1 481898.8 535820.2
## 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.1
## 2025 507682.3 547763.4 592944.4 524602.7Variable 2
# Convertir/declarar el PNCAFE en serie de tiempo mensual
variable2_ts <- ts(data_col$PNCAFE, start = c(2012, 1), frequency = 12)
variable2_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
## 2025 1355.7672 1361.3928 1063.9384 702.8167 818.6413 909.1030 1373.3560
## 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.2306
## 2025 1242.7417 1142.4126 1208.0865 1266.0696 1233.4232Variable 3
# Convertir/declarar el PECAFE en serie de tiempo mensual
variable3_ts <- ts(data_col$PECAFE, start = c(2012, 1), frequency = 12)
variable3_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
## 2025 345.9457 408.1100 404.0800 390.2119 401.8219 366.9740 323.6595 352.6757
## 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.5167
## 2025 406.1490 401.6996 407.6458 384.9918CONTEXTO
Sector económico elegido:
*Agricultura - CAFE*
Empresa dentro del Sector: Juan Valdez
Es una empresa colombiana de cafeterías y productos de café creada por la Federación Nacional de Cafeteros de Colombia. Actualmente cuenta con alrededor de 633 de tiendas en Colombia y en varios países, y es una de las marcas más reconocidas del café colombiano en el mundo. Su objetivo es promover y comercializar el café 100 % colombiano producido por miles de familias caficultoras del país.
Justificación de variables:
1.La Producción Nacional de Café (PNCAFE): Representa la oferta agregada del sector y determina la disponibilidad del grano tanto para el mercado interno como para la exportación. Un aumento en la producción suele traducirse en una mayor capacidad exportadora y abastecimiento interno, lo que puede estimular el crecimiento del sector.
3.Precio Externo del Café Colombiano (PECAFE): Captura la cotización del café en los mercados internacionales. Debido a la fuerte integración del café colombiano al comercio global, las variaciones en este precio influyen directamente en los ingresos por exportaciones y en los incentivos para la producción.
4.Exportaciones de Café (XCAF): Permite evaluar el desempeño del sector en el comercio internacional y reflejan la demanda externa por el café colombiano. Esta variable está estrechamente vinculada tanto a la producción nacional como al comportamiento de los precios internacionales.
METODOLOGÍA
Para analizar el comportamiento de las series temporales utilizadas en el estudio, se empleó la técnica de descomposición STL (Seasonal and Trend decomposition using Loess). Este método permite separar una serie temporal en tres componentes principales: tendencia, estacionalidad y residuo.
La tendencia refleja el comportamiento de largo plazo de la serie, mostrando la dirección general de la variable a lo largo del tiempo. La estacionalidad captura los patrones que se repiten de manera periódica en intervalos regulares, como variaciones mensuales o anuales. Además, el componente residual o irregular recoge las fluctuaciones aleatorias o choques inesperados que no pueden ser explicados por la tendencia ni por la estacionalidad. Tambien, el método STL se basa en el procedimiento de suavizamiento, el cual permite estimar de manera flexible los componentes de la serie sin imponer supuestos estrictos sobre su forma funcional. Debido a esta flexibilidad, la descomposición STL es ampliamente utilizada en el análisis de series temporales económicas y agrícolas.
Asimismo, se aplicó el modelo ARIMA siguiendo la metodología Box–Jenkins. Este enfoque se desarrolla en cuatro etapas: identificación, estimación, validación y pronóstico. En la primera fase se evaluó la estacionariedad de la serie mediante gráficos y la prueba de Dickey–Fuller, aplicando diferenciaciones cuando fue necesario. Posteriormente, se analizaron las funciones de autocorrelación y autocorrelación parcial para determinar los posibles valores de los parámetros autorregresivos y de medias móviles. Una vez definidos los órdenes 𝑝,𝑑,𝑞, se procedió a la estimación de los coeficientes mediante máxima verosimilitud, verificando su significancia estadística. La validación del modelo incluyó pruebas de Ljung–Box y Durbin Watson para confirmar que los residuos se comportaran como ruido blanco, además de la comparación de criterios de información como AIC y BIC. Finalmente, se evaluó la capacidad predictiva del modelo separando la muestra en entrenamiento y prueba, calculando métricas de error como RMSE y MAPE.
De manera adicional, se utilizó el procedimiento automático de selección de modelos ARIMA con el fin de corroborar los resultados obtenidos en la estimación manual. Este método permitió confirmar la validez del modelo elegido y reforzar la consistencia del análisis.
EXTRACCIÓN DE SEÑALES
La descomposición de series de tiempo permite comprender mejor los datos, mejorar predicciones y tomar decisiones más estratégicas. Es una herramienta clave en la analítica de negocios, especialmente en entornos donde las fluctuaciones en los datos pueden afectar inversiones, políticas económicas y estrategias empresariales.
Gráfico inicial de la variable 1 en niveles -Original
library(ggplot2)
library(plotly)
# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable1_ts)
# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)),
y = variable1)) +
geom_line(color = "grey", linewidth = 0.4) + # Cambiado 'size' por 'linewidth'
geom_point(color = "black", size = 0.1) +
ggtitle("Variable 1: Serie original") +
xlab("Tiempo") +
ylab("Unidad Variable 1") +
theme_minimal()
ggplotly(grafico_serie)Extracción señales variable 1 - XCAF:
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp_var1 <- stl(variable1_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var1 <- data.frame(
Time = rep(time(variable1_ts), 4), # Tiempo repetido para cada componente (son 4 componentes)
Value = c(stl_decomp_var1$time.series[, "seasonal"],
stl_decomp_var1$time.series[, "trend"],
stl_decomp_var1$time.series[, "remainder"],
variable1_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable1_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_var1, aes(x = Time, y = Value, color = Component)) +
geom_line() +
facet_wrap(~Component, scales = "free_y", ncol = 1) +
theme_minimal() +
labs(title = "Descomposición temporal de la variable 1",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Extracción señales variable 2 - PNCAFE
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp_var2 <- stl(variable2_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var2 <- data.frame(
Time = rep(time(variable2_ts), 4), # Tiempo repetido para cada componente
Value = c(stl_decomp_var2$time.series[, "seasonal"],
stl_decomp_var2$time.series[, "trend"],
stl_decomp_var2$time.series[, "remainder"],
variable2_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable2_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_var2, aes(x = Time, y = Value, color = Component)) +
geom_line() +
facet_wrap(~Component, scales = "free_y", ncol = 1) +
theme_minimal() +
labs(title = "Descomposición temporal de la variable 2",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Extracción señales variable 3 - PECAFE
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp_var3 <- stl(variable3_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_var3 <- data.frame(
Time = rep(time(variable3_ts), 4), # Tiempo repetido para cada componente
Value = c(stl_decomp_var3$time.series[, "seasonal"],
stl_decomp_var3$time.series[, "trend"],
stl_decomp_var3$time.series[, "remainder"],
variable3_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(variable3_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, aes(x = Time, y = Value, color = Component)) +
geom_line() +
facet_wrap(~Component, scales = "free_y", ncol = 1) +
theme_minimal() +
labs(title = "Descomposición temporal de la variable 3",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Después de la descomposición temporal de cada variable, se extrae la variable ajustada por estacionalidad para graficarla junto con la serie original:
Se crea la variable1 ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]Se crea la variable2 ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]Se crea la variable3 ajustada por estacionalidad
# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]Ahora si se puede graficar las series originales versus la ajustada por estacionalidad
Serie Original vs Ajustada por estacionalidad
Gráfico serie original VS ajustada Variable 1 - XCAF
# Crear vector de fechas correctamente alineado con la serie
fechas_var1 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var1 <- ggplot() +
geom_line(aes(x = fechas_var1, y = variable1_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Variable 1: Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Unidad de medida variable 1") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)Gráfico serie original VS ajustada Variable 2 - PNCAFE
# Crear vector de fechas correctamente alineado con la serie
fechas_var2 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var2 <- ggplot() +
geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Variable 2: Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Unidad de medida variable 2") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)Gráfico serie original VS ajustada Variable 3 - PECAFE
# Crear vector de fechas correctamente alineado con la serie
fechas_var3 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var3 <- ggplot() +
geom_line(aes(x = fechas_var3, y = variable3_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Variable 3: Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Unidad de medida variable 3") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)Serie original vs tendencia
La extracción de la tendencia permite centrarse en los cambios estructurales de la serie.
Analizar la tendencia ayuda a prever escenarios futuros y anticipar posibles crisis o oportunidades en el sector o variable de análisis
Primero se debe obtener la tendencia de cada variable y luego graficarla
Tendencia Variable 1 - XCAF
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Variable 1: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Unidad de medida Variable 1") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var1)Interpretación:
En este caso la tendencia muestra un cambio estructural: tras un periodo de relativa estabilidad hasta 2019, a partir de 2020 se observa un salto significativo en el nivel de exportaciones, que se mantiene en ascenso hasta 2024, a pesar de que hay una caída para principios del año 2025, logra recuperaste de manera positiva y acelerada manteniendo una tendencia positiva posiblemente hacia los próximos años. Esto no solo indica crecimiento, la tendencia refleja una transformación en la dinámica exportadora, pasando de un patrón estable a uno de expansión acelerada, lo que sugiere un reposicionamiento estratégico de Juan Valdez en el comercio global del café.
Tendencia Variable 2 - PNCAFE
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Variable 2: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Unidad de medida Variable 2") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var2)Interpretación:
La tendencia de la variable producción de Café, evidencia un crecimiento sostenido desde 2013, con una desaceleración entre los años 2021 hasta 2023 aproximadamente, pero se logra recuperar en el 2024, este comportamiento nos orienta a pensar que en los próximos años la empresa podría mantener una tendencia positiva si logra consolidar estrategias para los factores coyunturales.
Tendencia Variable 3 - PECAFE
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Variable 3: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Unidad de medida Variable 3") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var3)Interpretación:
El precio externo del café colombiano presenta un comportamiento muy similar al de las exportaciones de café: ambas variables muestran una relativa estabilidad para luego tener un crecimiento acelerado que se podría prolongar hacia el 2026. La tendencia revela una clara revalorización del café colombiano en el exterior en los últimos años.
Tasa de crecimiento anual de las variables
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia:
Tasa de crecimiento de la serie de tendencia y original para la variable 1 - XCAF
#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] 156print(length(tasa_crecimiento_var1))## [1] 156print(length(tasa_tendencia_var1))## [1] 156*Gráfico variable original y tendencia variable 1: tasa de crecimiento anual**
library(ggplot2)
library(plotly)
# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable1: Tasa de crecimiento anual % de la serie Original y la tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)Interpretación:
En este caso, la tasa de crecimiento revela un patrón más claro: después de un periodo de relativa estabilidad, a partir de 2020 se observa un crecimiento sostenido en la tasa de expansión, lo que indica un cambio en la dinámica exportadora. Este comportamiento es consistente con el aumento de los precios externos y refleja que Colombia, y en particular empresas como Juan Valdez, han logrado consolidar su presencia en mercados internacionales. Finalmente, aunque la tasa de crecimiento es un poco volátil en el corto plazo, la tendencia apunta hacia una expansión firme y estratégica de las exportaciones en el corto plazo.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2 - PNCAFE
#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] 156print(length(tasa_crecimiento_var2))## [1] 156print(length(tasa_tendencia_var2))## [1] 156# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable2: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)Interpretación:
La tasa de crecimiento anual de la producción de café colombiano evidencia un comportamiento muy volátil en la serie original, con picos de expansión y caídas abruptas que reflejan la sensibilidad del sector a factores coyunturales como el clima, los costos de insumos y la dinámica de la demanda externa. Sin embargo, la tendencia suavizada muestra un patrón de desaceleración progresiva hacia el final del periodo analizado, lo que se esperaría el mismo comportamiento para los siguientes años.
Este descenso en la tasa de crecimiento no significa necesariamente una caída en la producción, sino una moderación del ritmo de expansión. En otras palabras, la producción sigue aumentando, pero cada vez a un ritmo menor. Esto sugiere que el sector podría estar entrando en una fase de maduración, donde la estabilidad y la consolidación pesan más que el crecimiento acelerado, lo que obliga a replantear estrategias de productividad y competitividad para sostener la dinámica en el largo plazo.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 3 - PECAFE
#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] 156print(length(tasa_crecimiento_var3))## [1] 156print(length(tasa_tendencia_var3))## [1] 156# Gráfico de la tasa de crecimiento anual variable 3
grafico_crecimiento_var3 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable3: Tasa de crecimiento anual % de la serie Original y la tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)Interpretación:
La tasa de crecimiento anual del precio externo del café colombiano refleja una dinámica marcada por la volatilidad propia de los mercados internacionales: la serie original muestra picos de expansión en años como 2015, 2020 y 2023, intercalados con caídas que evidencian la sensibilidad del precio frente a factores globales como oferta mundial, especulación financiera y variaciones en la demanda.
La tendencia suavizada permite ver con mayor claridad que, pese a esas oscilaciones, el precio externo ha mantenido un comportamiento de crecimiento en el largo plazo. Este patrón es consistente con lo observado en las exportaciones de café, ya que ambos indicadores muestran un cambio estructural a partir de 2020: mientras las exportaciones aumentan en volumen, el precio externo se revaloriza, lo que sugiere un entorno internacional a corto plazo favorable para la competitividad del café colombiano y para empresas como Juan Valdez.
MODELO ARIMA
División en conjunto de entrenamiento y prueba para la variable 1 que es la elegida para pronosticar:
# Esta división idealmente podria se 80%-70% de los datos para entrenamiento y 20%-30% para prueba o test
# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 2025 y el conjunto de prueba o test: noviembre 2025-diciembre 2025
train_size <- length(variable1_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable1_ts, end = c(2025, 9)) # Entrenamiento hasta septiembre 2025
test_ts <- window(variable1_ts, start = c(2025, 10)) # Prueba inicia desde oct2025Paso 1: Identificación del modelo
Identificar estacionariedad
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
A continuación se aplica el test ADF para validar estacionariedad en el conjunto de entrenamiento de la variable 1, que es la elegida para pronosticar:
library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable 1 (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test##
## Augmented Dickey-Fuller Test
##
## data: train_ts
## Dickey-Fuller = -1.9592, Lag order = 5, p-value = 0.5933
## alternative hypothesis: stationaryEl test ADF en la variable 1 arrojó un p-value igual a 0.5933, este valor es mayor a 0.05, por tanto la serie no es estacionaria. Asimismo, esto se confirma al visualizar la gráfica de la serie original. De ese modo se debe ejecutar el código siguiente para diferenciar una vez la variable 1 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 variable 1
A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable 1 - PNCAFE para ver graficamente el cambio o ajuste:
# Graficar la serie original:
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "blue") +
ggtitle("Variable 1: Serie Original") +
xlab("Tiempo") + ylab("Gwh")
ggplotly(p2) # Convertir en gráfico interactivo # Graficar la serie diferenciada:
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("variable 1: Serie Estacionaria (Una diferenciación)") +
xlab("Tiempo") + ylab("variable 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoA continuación se aplica la diferenciación logarítimica y la varible XCAF u objeto ahora se llama train_diff_log:
# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación
train_diff_log <- diff(log(train_ts), differences = 1)La aplicación de logaritmos en series de tiempo se realiza principalmente para lograr que la serie sea más estable, lineal y estacionaria. Esta transformación es relevante porque permite modelar mejor las series que siguen un crecimiento exponencial y facilita la aplicación de técnicas estadísticas que requieren estacionariedad.
Ahora graficamos la serie orignal versus la serie diferenciada una vez con logaritmo
# Graficar la serie original:
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
geom_line(color = "blue") +
ggtitle("Variable1: Serie Original") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p2) # Convertir en gráfico interactivo # Graficar la serie diferenciada en logaritmo:
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], micro_Diff = as.numeric(train_diff_log)), aes(x = Tiempo, y = micro_Diff)) +
geom_line(color = "red") +
ggtitle("Variable1: Serie Estacionaria (Una diferenciación en logaritmo)") +
xlab("Tiempo") + ylab("Variable 1 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoAhora probamos estacionariedad en la serie diferenciada ( en nivel y logaritmo)
# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
adf_test_diff <- adf.test(train_diff)
print(adf_test_diff)##
## Augmented Dickey-Fuller Test
##
## data: train_diff
## Dickey-Fuller = -7.159, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary# Segunda prueba de estacionariedad sobre la serie diferenciada en logaritmo
adf_test_diff_log <- adf.test(train_diff_log)
print(adf_test_diff_log)##
## Augmented Dickey-Fuller Test
##
## data: train_diff_log
## Dickey-Fuller = -7.6693, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEn el test ADF se muestra que el valor que puede tomar d=1:
El p-value ya es menor a 0.05 con una primera diferencia en ambos casos: niveles o con logaritmo natural. Por tanto el valor que puede tomar d es igual a 1.
Identificación manual de p y q
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_log, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")
ggplotly(acf_plot)ggplotly(pacf_plot)Interpretación correlogramas
Se puede observar que los valores que podrian tomar p y q serian:
p=1,2,3,6 y q=1 (P optimo=6) (q óptimo=1)
Posibles Modelos (p,d,q) (1,1,1)
(2,1,1)
(3,1,1) (6,1,1)
El modelo óptimo para esta variable seria (6,1,1)
Paso 2. Estimación manual del modelo
Estimación del modelo identificado (6,1,1)
# Cálculo manual de modelo ARIMA
manual_arima_model_1 <- Arima(train_ts, order = c(6,1,1)) #Se va a estimar un modelo Arima de orden (6,1,1)
summary(manual_arima_model_1)## Series: train_ts
## ARIMA(6,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ma1
## -0.3341 -0.2974 -0.1086 -0.1080 -0.0710 -0.1212 -0.1212
## s.e. 0.3662 0.1806 0.1496 0.0991 0.0888 0.0797 0.3635
##
## sigma^2 = 2.356e+09: log likelihood = -1998.9
## AIC=4013.81 AICc=4014.74 BIC=4038.61
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4301.009 47343.97 35036.9 -1.818374 15.5157 0.5465585 -0.0123449Significancia de coefientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model_1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.334102 0.366185 -0.9124 0.36157
## ar2 -0.297389 0.180620 -1.6465 0.09966 .
## ar3 -0.108609 0.149621 -0.7259 0.46790
## ar4 -0.108030 0.099074 -1.0904 0.27554
## ar5 -0.071038 0.088829 -0.7997 0.42387
## ar6 -0.121180 0.079710 -1.5203 0.12845
## ma1 -0.121165 0.363513 -0.3333 0.73890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación significancia coeficientes
Los coeficiente parecen no tener un impacto importante en el modelo, debido a que no son significativos al 5%. Sin embargo, con un nivel de significancia al 10%, podemos identificar que almenos uno de los componentes (ar2) es 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.
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_1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(6,1,1)
## Q* = 33.807, df = 17, p-value = 0.00889
##
## Model df: 7. Total lags used: 24Serie de residuos (gráfico superior)
Muestra cómo se comportan los errores a lo largo del tiempo. Se observan picos, lo que podría indicar cambios estructurales en los datos.
ACF de residuos La mayoría de las barras están dentro de las líneas azules (intervalos de confianza). Sin embargo, algunos rezagos parecen salir del rango, lo que sugiere que aún puede haber estructura no capturada en los datos. Esto indica que el modelo podría mejorarse pero es un modelo aceptable.
Histograma de residuos con ajuste normal Los residuos se acercan a la normalidad, pero hay algunos valores extremos (colas más gruesas de lo esperado). Esto indica que puede haber eventos atípicos o datos no bien explicados por el modelo.
Paso 4. Pronóstico (modelo manual)
Pronóstico en el test de prueba (oct, nov y dic 2025) 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_1, 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("Variable1: Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p_manual)El modelo (6,1,1) parece tener un sub ajuste, pero se capturan bien los puntos de quiebre. Es un modelo tentativo adecuado para pronostico fuera de muestra o a futuro.
Métricas de evaluación del modelo manual dentro del periodo de prueba (oct,nov y dic2025
# 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: 49725.73cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 55931.59A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic2025
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_manual <- forecast(manual_arima_model_1, 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 2025.750 547763.4 504888.8
## 2 2025.833 592944.4 508998.0
## 3 2025.917 524602.7 502246.5Ahora pronosticamos fuera del periodo de análisis: Enero 2026
Es decir, le sumamos al periodo de prueba (oct,nov,dic2025) una observación más (enero 2026). 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_1, 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 2025.750 504888.8
## 2 2025.833 508998.0
## 3 2025.917 502246.5
## 4 2026.000 495276.6# 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: 2026 = 495276.559509132"El modelo (6,1,1) logra captar el quiebre y se podría considerar que representa un pronóstico válido, aunque con valores pronósticados subestimados.
Estimación de los otros modelos posibles
Adicionalmente, probamos con los demás posibles modelos: (1,1,1), (2,1,1) y (3,1,1)
Estimación del modelo (1,1,1)
# Cálculo manual de modelo ARIMA
manual_arima_model_2 <- Arima(train_ts, order = c(1,1,1)) #Se va a estimar un modelo Arima de orden (1,1,1)
summary(manual_arima_model_2)## Series: train_ts
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.1755 -0.6435
## s.e. 0.1343 0.1016
##
## sigma^2 = 2.328e+09: log likelihood = -2000.44
## AIC=4006.89 AICc=4007.04 BIC=4016.19
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 4214.328 47803.83 35331 -1.884619 15.60854 0.5511463 -0.002035483Significancia de coefientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model_2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.17550 0.13433 1.3065 0.1914
## ma1 -0.64347 0.10164 -6.3312 2.433e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 checkresiduals(manual_arima_model_2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)
## Q* = 35.161, df = 22, p-value = 0.03727
##
## Model df: 2. Total lags used: 24Pronóstico en el test de prueba (oct, nov y dic 2025) 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_2, 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("Variable1: Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p_manual)Estimación del modelo (2,1,1)
# Cálculo manual de modelo ARIMA
manual_arima_model_3 <- Arima(train_ts, order = c(2,1,1)) #Se va a estimar un modelo Arima de orden (2,1,1)
summary(manual_arima_model_3)## Series: train_ts
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 0.0506 -0.1154 -0.5097
## s.e. 0.2235 0.1257 0.2189
##
## sigma^2 = 2.329e+09: log likelihood = -2000.01
## AIC=4008.02 AICc=4008.27 BIC=4020.42
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3988.317 47675.95 35342.04 -1.946048 15.68955 0.5513185
## ACF1
## Training set -0.006283686Significancia de coefientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model_3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.05063 0.22351 0.2265 0.82080
## ar2 -0.11539 0.12574 -0.9177 0.35877
## ma1 -0.50970 0.21887 -2.3287 0.01987 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 checkresiduals(manual_arima_model_3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1)
## Q* = 33.909, df = 21, p-value = 0.03706
##
## Model df: 3. Total lags used: 24Pronóstico en el test de prueba (oct, nov y dic 2025) 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_3, 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("Variable1: Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p_manual)Estimación del modelo (3,1,1)
# Cálculo manual de modelo ARIMA
manual_arima_model_4 <- Arima(train_ts, order = c(3,1,1)) #Se va a estimar un modelo Arima de orden (3,1,1)
summary(manual_arima_model_4)## Series: train_ts
## ARIMA(3,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ma1
## 0.1823 -0.0604 0.0688 -0.6389
## s.e. 0.2705 0.1404 0.1248 0.2561
##
## sigma^2 = 2.34e+09: log likelihood = -1999.86
## AIC=4009.73 AICc=4010.11 BIC=4025.23
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 4101.568 47632.59 35278.27 -1.917475 15.6628 0.5503236
## ACF1
## Training set -0.008725256Significancia de coefientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model_4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.182307 0.270542 0.6739 0.50040
## ar2 -0.060401 0.140371 -0.4303 0.66698
## ar3 0.068815 0.124831 0.5513 0.58145
## ma1 -0.638911 0.256104 -2.4947 0.01261 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 checkresiduals(manual_arima_model_4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,1,1)
## Q* = 34.502, df = 20, p-value = 0.02292
##
## Model df: 4. Total lags used: 24Pronóstico en el test de prueba (oct, nov y dic 2025) 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_4, 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("Variable1: Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable1")
ggplotly(p_manual)Al aplicar forecast, podemos confirmar que los modelos (1,1,1), (2,1,1) y (3,1,1) no logran captar de manera adecuara los puntos de quiebre y tampoco logran pronósticos válidos
MODELO ARIMA AUTOMÁTICO
Identificación automática del modelo ARIMA
library(forecast)
# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)
# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)## Series: train_ts
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.8485 0.3959 -0.5332
## s.e. 0.0628 0.0806 0.0683
##
## sigma^2 = 2.275e+09: log likelihood = -1998.25
## AIC=4004.5 AICc=4004.75 BIC=4016.9
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3932.893 47119.03 33776.97 -1.850663 14.93934 0.5269042 0.01120106Significancia de coeficientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.848507 0.062776 -13.5164 < 2.2e-16 ***
## ma1 0.395889 0.080598 4.9119 9.02e-07 ***
## ma2 -0.533167 0.068283 -7.8082 5.80e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación significancia coeficientes: ar1, ma1, ma2 son altamente significativos, lo que indica que todos los coeficientes tienen un impacto importante en el modelo.
# Ajuste del modelo ARIMA(1,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(1, 1, 2))
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.8485 0.3959 -0.5332
## s.e. 0.0628 0.0806 0.0683
##
## sigma^2 = 2.275e+09: log likelihood = -1998.25
## AIC=4004.5 AICc=4004.75 BIC=4016.9
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3932.893 47119.03 33776.97 -1.850663 14.93934 0.5269042 0.01120106El modelo automático identificado es (1,1,2). Si se compara el AIC o BIC de este modelo frente el modelo manual (6,1,1), se obtiene un valor más bajo en esta métricas en este modelo automático. Probablemente pudiera ser un buen modelo para pronosticar la variable 1 - XCAF.
# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_auto) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2)
## Q* = 34.827, df = 21, p-value = 0.0295
##
## Model df: 3. Total lags used: 24Serie de residuos (gráfico superior)
Se observan picos, lo que podría indicar cambios estructurales en los datos.
ACF de residuos La mayoría de las barras están dentro de las líneas azules. Asimismo, se observa que solo hay pocos rezagos parecen salir del rango, lo que que aún puede haber una menor parte de estructura no capturada en los datos (a comparación del modelo (6,1,1)).
Histograma de residuos con ajuste normal Los residuos se acercan a la normalidad, hay algunos valores extremos. Esto indica que puede haber algunos eventos atípicos o datos no bien explicados por el modelo.
Pronóstico modelo ARIMA automático (1,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 interactivoEn el modelo automático (1,1,2) hay un sub ajuste, pero se capturan muy bien los puntos de quiebre. Es un modelo tentativo adecuado para pronpostico fuera de muestra o a futuro.
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table_auto <- data.frame(
Tiempo = time(arima_forecast_auto$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast_auto$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table_auto)## Tiempo Observado Pronosticado
## 1 2025.750 547763.4 493522.0
## 2 2025.833 592944.4 501971.2
## 3 2025.917 524602.7 494802.0Ahora 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 2025.750 493522.0
## 2 2025.833 501971.2
## 3 2025.917 494802.0
## 4 2026.000 500885.1# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2026:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2026: 2026 = 500885.115420756"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar el mes enero 2026
future_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = 4)
# Extraer el valor específico de enero 2026
forecast_ener2026_auto <- future_forecast_auto$mean[4]
print(paste("Pronóstico para enero 2026:", forecast_ener2026_auto))## [1] "Pronóstico para enero 2026: 500885.115420756"MODELO SARIMA AUTOMÁTICO
El modelo SARIMA(1,1,2)(0,0,1)[12] sugiere que:
- La serie tiene una tendencia no estacionaria, corregida con una diferenciación.
- Hay un componente estacional autorregresivo fuerte cada 12 períodos.
- El ajuste es adecuado según los criterios AIC y BIC, pero se podría comparar con otros modelos para mejorar la predicción.
# Identificación automática del 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)(0,0,1)[12]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## -0.8728 0.4146 -0.5324 0.1545
## s.e. 0.0544 0.0762 0.0676 0.0753
##
## sigma^2 = 2.23e+09: log likelihood = -1996.18
## AIC=4002.37 AICc=4002.74 BIC=4017.86A 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)(0,0,1)[12]
darima <- Arima(train_ts,
order = c(1, 1, 2), # (p,d,q) -> (0,1,1)
seasonal = list(order = c(0, 0, 1), # (P,D,Q) -> (1,0,0)
period = 12)) # Periodicidad estacional de 12 meses
# Mostrar resumen del modelo ajustado
summary(darima)## Series: train_ts
## ARIMA(1,1,2)(0,0,1)[12]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## -0.8728 0.4146 -0.5324 0.1545
## s.e. 0.0544 0.0762 0.0676 0.0753
##
## sigma^2 = 2.23e+09: log likelihood = -1996.18
## AIC=4002.37 AICc=4002.74 BIC=4017.86
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3705.497 46498.56 33950.37 -1.818631 15.0316 0.5296091 0.001637645En el correlograma de residuos siguiente se observa que, mejora la correlación de los residuos frente a los dos modelos anteriores. Sin embargo, el modelo no captura de manera adecuada los datos para realizar pronósticos
# 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)(0,0,1)[12]
## Q* = 27.377, df = 20, p-value = 0.125
##
## Model df: 4. 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 interactivoAl aplicar forecast, podemos confirmar que el modelo SARIMA no logra captar de manera adecuara el punto de quiebre y tampoco logran pronósticos válidos. (SARIMA NO ES UN MODELO VÁLIDO EN ESTE CASO)
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(auto_arima_model, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
Tiempo = time(arima_forecast$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table)## Tiempo Observado Pronosticado
## 1 2025.750 547763.4 502617.4
## 2 2025.833 592944.4 506938.1
## 3 2025.917 524602.7 520137.1Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueba una observación más. Es decir, se estan pronosticando 4 observaciones o meses.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
Tiempo = time(next_forecast$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast)## Tiempo Pronostico
## 1 2025.750 502617.4
## 2 2025.833 506938.1
## 3 2025.917 520137.1
## 4 2026.000 529353.7# 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: 2026 = 529353.664622314"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)(0,0,1)[12]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## -0.8728 0.4146 -0.5324 0.1545
## s.e. 0.0544 0.0762 0.0676 0.0753
##
## sigma^2 = 2.23e+09: log likelihood = -1996.18
## AIC=4002.37 AICc=4002.74 BIC=4017.86# 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: 502617.395208676"CONCLUSIÓN
El modelo automático ARIMA(1,1,2) fue seleccionado como el más adecuado para el pronóstico de la variable XCAF, al presentar un mejor desempeño relativo frente al modelo manual ARIMA(6,1,1). En particular, este modelo registra valores más bajos en los criterios de información AIC y BIC, lo que sugiere un mejor ajuste y mayor capacidad predictiva. Asimismo, sus coeficientes resultan estadísticamente significativos al 5%, evidenciando que las variables incluidas tienen un impacto relevante en la dinámica de la serie.
A partir de este modelo, se obtuvo un pronóstico para enero de 2026 de 500,885.1 miles de dólares, acompañado de una tasa de crecimiento anual de 6%, lo que indica que las exportaciones de café se mantendrían en terreno positivo, aunque con señales de desaceleración. No obstante, la tasa de crecimiento mensual de −5% refleja una contracción en el corto plazo, evidenciando una pérdida de dinamismo reciente en el comportamiento de la serie.
Si bien se identifica la presencia de cierto subajuste en el modelo automático, este logra capturar adecuadamente los puntos de quiebre, lo cual es clave para fines de pronóstico. No obstante, el análisis de los residuos sugiere que aún existen oportunidades de mejora en la especificación del modelo.
EVALUACIÓN DEL MODELO
Las exportaciones de café colombiano podrían mantener un crecimiento anual positivo, pero con una desaceleración, lo cuál se coincide con la tasa de crecimiento anual pronosticada (6%) para enero 2026. Esto se puede relacionar con:
La caída sostenida en la tasa de crecimiento de la PNCAFE, que limita la expansión del volumen exportado, restringiendo así el crecimiento de las exportaciones. Este menor dinamismo por el lado de la oferta refuerza la desaceleración observada en la tasa anual y puede explicar reducciones puntuales de los envíos en el corto plazo, coherentes con la variación mensual negativa para enero 2026.
El PECAFE sigue aumentando, pero cada vez menos (desaceleración), entonces el ingreso total por exportaciones también crece, pero a un ritmo menor. Aunque siga exportando y obteniendo más ingresos, el hecho de que los precios ya no suban con la misma fuerza hace que ese aumento sea más moderado, lo que se refleja en el 6% anual (positivo, pero desacelerado).
IMPLICACIONES DEL PRONÓSTICO
El pronóstico del 6% anual y el -5% mensual de las exportaciones de café, funcionan como un indicador del entorno externo en el que opera Juan Valdez.
En este sentido, tendría las siguientes implicaciones:
Expansión en el mercado internacional: El crecimiento anual positivo pero desacelerado (6%) se traduce en un entorno donde la demanda externa o los ingresos por exportaciones crecen más lentamente, lo que puede limitar el ritmo al que la empresa expande sus ventas internacionales.
Reducción de ingresos y ventas externas: La caída mensual de - 5% refleja un ajuste de corto plazo en las exportaciones nacionales. Esto afecta a la empresa en la medida en que forma parte de ese agregado de exportaciones, implicando menores ingresos por exportaciones en el corto plazo y/o mayor volatilidad en sus ventas externas.
RECOMENDACIONES ESTRATÉGICAS
Decisiones
La empresa debe fortalecer su capacidad productiva para sostener el crecimiento de las exportaciones y aprovechar la revalorización del precio externo. Esto implica invertir en productividad agrícola, innovación tecnológica y diversificación de mercados. Además, debe consolidar estrategias de diferenciación de calidad para mantener precios competitivos en el largo plazo.
Riesgos y oportunidades
Riesgos: La desaceleración en la tasa de crecimiento de la producción nacional limita la oferta exportable; la volatilidad de precios internacionales y la dependencia de factores externos (clima, costos logísticos) pueden afectar ingresos.
Oportunidades: La tendencia positiva en exportaciones y precios externos abre espacio para expandir presencia internacional, consolidar alianzas comerciales y posicionar el café colombiano como producto premium.
Estrategias
Diversificar mercados y canales de exportación para reducir riesgos de concentración.
Programas con responsabilidad social con apoyo especialmente a los agricultores.
Invertir en sostenibilidad y certificaciones que refuercen la competitividad internacional.
Se debe enfocarse en aumentar productividad, aprovechar la valorización externa, asegurando un crecimiento estable en el futuro cercano.