Instalar/Cargar librerias necesarias para el análisis

#Cargar librerías necesarias
library(readxl)  # Para leer archivos Excel
library(tseries)  # Para pruebas de estacionariedad
library(forecast)  # Para modelado ARIMA y pronósticos
library(ggplot2)  # Para visualización de datos
library(plotly)  # Para gráficos interactivos
library(timetk)   #timetk simplifica y acelera el análisis exploratorio, visualización, y preparación de datos temporales para modelado. Es ideal para quienes trabajan con series temporales en un flujo de trabajo "tidy" y buscan integrar análisis visuales, detección de patrones y forecasting en un solo paquete.

Cargar base de datos

library(readxl)
data_col <- read_excel("~/Desktop/Analisis de Negocios/Caso 2. CAFE/Base Caso2 (2).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"))

Introducción

¿Por qué elegiste este sector y empresa?

Elegimos este sector ya que el cafe es uno de los productos más representativos de Colombia, con alto impacto en la economía nacional y en la generación de divisas.

La empresa exportadora de café se ubica en un sector estratégico, pues conecta la producción agrícola con los mercados internacionales, lo que la hace relevante para inversionistas y directivos.

Además, es un sector con datos históricos amplios y confiables, lo que permite aplicar modelos estadísticos robustos

Por otro lado, elegimos la empresa ficticia Cafexport SAS con cultivos en en el Cauca y Eje Cafetero, cuyos mercados principales son:

Estados Unidos: principal destino, con contratos de largo plazo.

Europa: Bélgica, Alemania y Francia, enfocados en cafés premium.

Asia: Japón y Corea del Sur, interesados en cafés sostenibles y de origen único.

Y los retos actuales que enfrentan son:

Volatilidad climática: fenómenos de El Niño y La Niña afectan la producción.

Costos logísticos: transporte internacional y disponibilidad de contenedores.

Competencia global: Brasil y Vietnam como grandes productores.

Seguridad rural: presiones sociales en zonas cafeteras.

¿Por qué son relevantes las tres variables?

Producción de café: mide la capacidad de oferta y refleja la resiliencia del sector frente a factores climáticos y sociales.

Precio externo del café colombiano : determina los ingresos potenciales en divisas y la competitividad frente a otros países productores.

Exportaciones de café: muestran el desempeño comercial y la capacidad para posicionarse en mercados internacionales.

estas tres variables permiten construir un mapa integral del negocio cafetero, conectando la oferta agrícola con la dinámica internacional y ofreciendo insumos clave para decisiones estratégicas de producción, comercialización y gestión financiera.

¿Qué relación esperamos encontrar entre ellas?

Una relación directa entre producción y exportaciones: mayor producción debería permitir mayores volúmenes exportados.

Una relación inversa entre producción global y precio externo: cuando la oferta mundial cae (ej. por clima en Brasil), el precio externo sube, beneficiando a Colombia.

Una relación positiva entre precio externo y exportaciones: precios altos incentivan mayores ventas externas y contratos de largo plazo.

¿Qué decisión empresarial podría apoyarse en esta información? Planificación de producción: ajustar estrategias agrícolas según pronósticos de precios y demanda externa.

Gestión de riesgos: usar coberturas cambiarias y contratos futuros para mitigar la volatilidad del precio internacional.

Diversificación de mercados: identificar destinos estratégicos cuando la producción y precios favorecen la exportación.

Inversión en sostenibilidad: priorizar variedades resistentes y cafés especiales para aprovechar precios premium en el mercado externo.

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

Variable 1

# Convertir/declarar PNCAFE (Producción de Café) en serie de tiempo mensual
variable1_ts <- ts(data_col$PNCAFE, start = c(2012, 1), frequency = 12)

Variable 2

# Convertir/declarar PECAFE (Precio externo del Café) en serie de tiempo  mensual
variable2_ts <- ts(data_col$PECAFE, start = c(2012, 1), frequency = 12)

Variable 3

# Convertir/declarar XCAF (Exportaciones de cafe colombiano) en serie de tiempo mensual
variable3_ts <- ts(data_col$XCAF, start = c(2012, 1), frequency = 12)

Extracción de señales

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

library(ggplot2)
library(plotly)

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

# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)), 
                                      y = variable1)) +
  geom_line(color = "#D2B48C", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "#8B4513", size = 0.3) +
  ggtitle("Producción de Cafe: Serie original") +
  xlab("Tiempo") +
  ylab("Miles de sacos de 60 Kg de café verde equivalente") +
  theme_minimal()

ggplotly(grafico_serie)

Extracción señales variable 1

# 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 Producción de Café",
       x = "Tiempo",
       y = "Valor")

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

De acuerdo con la informacion del componente estacional para la variable de Produccion de Cafe (PNCAFE):

Podemos decir que el inicio de año siempre presenta valores negativos, explicado porque en Enero y Febrero los cafetales se encuentran en fase de desarrollo y maduración del fruto. En este periodo se realizan labores de mantenimiento agronómico como fertilización, control de plagas, poda y preparación de los cafetales para la cosecha secundaria. Aquí también influyen los factores climáticos, ya que a comienzos de año suelen presentarse lluvias o sequías que afectan el ciclo productivo. Los productores aprovechan este tiempo para planificar la logística de recolección y contratar la mano de obra necesaria para la mitaca de abril–junio.

Por el contrario, el pico más alto se observa en los meses de agosto y septiembre, explicado porque en este periodo se inicia la cosecha principal. El clima es determinante: lluvias moderadas favorecen la maduración, mientras que excesos pueden retrasar la recolección. La logística rural se activa al máximo: transporte, beneficio húmedo y secado trabajan con plena capacidad. Es el momento en que el país asegura gran parte de su oferta anual de café, consolidando su papel en el mercado internacional.

Extracción señales variable 2

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

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

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

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var2, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 2",
       x = "Tiempo",
       y = "Valor")

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

De acuerdo con la informacion del componente estacional para la variable de Precio Externo de Cafe (PECAFE):

En enero y febrero, los cafetales colombianos están en fase de desarrollo y maduración del fruto, no en cosecha. Esto significa que la oferta exportable es limitada, y la menor disponibilidad de café en el mercado internacional tiende a elevar los precios. Además:

  • Los compradores internacionales (tostadores y distribuidores) suelen reponer inventarios a comienzos de año para garantizar suministro.
  • La demanda se mantiene estable o incluso aumenta, mientras la oferta colombiana es baja, generando presión alcista.
  • Los costos logísticos internacionales (transporte marítimo, contenedores) suelen ser más altos en el primer trimestre, lo que se refleja en el precio final.

En cambio en los meses de Mayo a Julio se observa un descenso en los precios externos, las razones principales son: - Varios países productores (Brasil, Vietnam, Honduras) ya han iniciado o están próximos a sus cosechas principales, lo que aumenta la disponibilidad de café en el mercado internacional. - La mayor oferta presiona los precios a la baja, especialmente en el segmento de cafés estándar. - Mayo y julio no son meses de alta demanda internacional: los compradores ya han repuesto inventarios a comienzos de año y esperan la entrada de la cosecha principal de Colombia en septiembre–diciembre.

Lo anterior, genera un periodo de menos presion de compra, contribuyendo a precios mas bajos.

Extracción señales variable 3

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

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

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

# Crear gráfico con ggplot2
p <- ggplot(stl_df_var3, aes(x = Time, y = Value, color = Component)) +
  geom_line() +
  facet_wrap(~Component, scales = "free_y", ncol = 1) + 
  theme_minimal() +
  labs(title = "Descomposición temporal de la variable 3",
       x = "Tiempo",
       y = "Valor")

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

De acuerdo con la informacion del componente estacional para la variable de Exportaciones de Cafe (XCAF):

Los menores niveles de exportación se registran en abril, mes que coincide con el inicio de la cosecha secundaria (mitaca), caracterizada por un volumen mucho menor que la cosecha principal de septiembre–diciembre. En este periodo, los cafetales apenas entran en fase de recolección, lo que limita la oferta exportable y se traduce en menores envíos al exterior. La demanda externa se mantiene estable, pero no ejerce presión alcista, ya que los importadores esperan la gran cosecha colombiana del segundo semestre para abastecerse.

En septiembre, las exportaciones de café colombiano suelen alcanzar sus niveles más altos, ya que coincide con el inicio de la cosecha principal. Los cafetales entran en fase de maduración masiva, generando un volumen significativo de café listo para recolección y esto se traduce en una mayor disponibilidad de café exportable, lo que impulsa los envíos al exterior. La logística rural y portuaria tambien se activa al máximo: transporte, beneficio húmedo y secado trabajan con plena capacidad y los exportadores aprovechan este periodo para cumplir contratos internacionales y sacar ventaja de la ventana de alta producción.

En conclusion, las tres variables muestran un patrón cíclico y altamente interdependiente:

  • El precio externo responde a la oferta y demanda global.
  • Las exportaciones dependen de la estacionalidad de la cosecha y la logística.
  • La producción nacional es el eje que condiciona ambas, afectada por clima y manejo agronómico.

El análisis evidencia que el sector cafetero colombiano está marcado por la volatilidad internacional y la estacionalidad interna, lo que exige estrategias de:

Diversificación de mercados, Cobertura cambiaria y logística eficiente, asi como Inversión en resiliencia agronómica.

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

Gráfico serie original VS ajustada Variable 1

# 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 = "#D2B48C", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var1, y = variable1_sa), color = "#8B4513", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("PNCAFE:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Miles de sacos de 60 Kg de café verde equivalente") +
  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

# 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 = "#D2B48C", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var2, y = variable2_sa), color = "#8B4513", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("PECAFE:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Centavos de dólar por libra de 453.6 gr") +
  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

# 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 = "#D2B48C", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var3, y = variable3_sa), color = "#8B4513", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("XCAF:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  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)

Ahora graficamos serie original vs tendencia

Primero se debe obtener la tendencia de cada variable y luego graficarla

Tendencia Variable 1

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
  geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "#D2B48C", "Tendencia" = "#8B4513")) +
  ggtitle("PNCAFE: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Miles de sacos de 60 Kg de café verde equivalente") +
  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)

Tendencia Variable 2

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
  geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "#D2B48C", "Tendencia" = "#8B4513")) +
  ggtitle("PECAF: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Centavos de dólar por libra de 453.6 gr") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

Tendencia Variable 3

library(ggplot2)
library(plotly)

# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])

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

# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
  geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
  geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
  scale_color_manual(values = c("Serie Original" = "#D2B48C", "Tendencia" = "#8B4513")) +
  ggtitle("XCAF: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  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)

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

#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] 156
print(length(tasa_crecimiento_var1))
## [1] 156
print(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 = "#D2B48C", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "#8B4513", size = 0.8, linetype = "dashed") +
  ggtitle("PNCAFE: 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)

De acuerdo a esta grafica, se evidencia que la variable no se mantiene estable, presentando fluctuaciones importantes a lo largo del tiempo.

Desde 2013 se observa un crecimiento positivo. Esto se explica porque Colombia venía saliendo de la crisis productiva causada por la roya (2008–2011). La Federación Nacional de Cafeteros promovió la renovación de cafetales con la variedad Castillo, resistente a la roya y con altos niveles de productividad. Gracias a esta estrategia, la producción nacional se recuperó rápidamente y se reflejó en el repunte inicial.

Entre agosto de 2013 y septiembre de 2014 se observa un descenso.La fase de La Niña trajo lluvias intensas en varias zonas cafeteras, afectando la floración y maduración del grano. Esto redujo la productividad y generó una disminución temporal en la oferta.Adicionalmente, los precios internacionales del café se mantuvieron bajos, lo que desincentivó la inversión en fertilización y mantenimiento

Desde finales de 2015 hasta inicios de 2018 la tasa de crecimiento desciende.El Fenómeno de El Niño (2015–2016) provocó sequías prolongadas que redujeron la floración y el llenado del grano. Posteriormente, La Niña (2017–2018) trajo lluvias excesivas, afectando la maduración y aumentando la incidencia de enfermedades como la roya. Esta combinación de sequías y lluvias intensas generó una caída sostenida en la productividad. En paralelo, la apreciación del dólar y la volatilidad de la TRM afectaron los ingresos de los productores, reduciendo su capacidad de reinversión.

En 2019 se observa una leve recuperación con crecimiento positivo, explicada por la maduración de los cafetales renovados y condiciones climáticas más estables. Sin embargo, en 2020 la tasa de crecimiento se desacelera y pasa a negativo, afectada por la pandemia, que generó problemas logísticos, escasez de mano de obra y retrasos en la recolección. Aunque la demanda internacional de café en el hogar aumentó, la producción nacional enfrentó serias dificultades operativas.

Entre 2021 y 2022, el Fenómeno de La Niña provocó lluvias prolongadas, reduciendo la productividad y generando picos irregulares. A esto se sumaron presiones de seguridad rural y la volatilidad de los precios internacionales, que alteraron los incentivos de producción y generaron incertidumbre en el sector.

En 2023 comienza un crecimiento sostenido, impulsado por condiciones climáticas más favorables, la consolidación de cafetales renovados y una mayor demanda internacional de cafés especiales y sostenibles. Este repunte se mantiene hasta buena parte de 2024.

Desde noviembre de 2024 la tasa de crecimiento vuelve a caer en negativo. La explicación está en una combinación de exceso de lluvias, menor inversión agronómica y mayores costos logísticos, que redujeron la capacidad productiva y exportadora.

Este comportamiento anticipa la fuerte contracción de 2025, sin embargo se considera un año mixto para el sector dado que si bien la producción tuvo una ligera disminución para el año 2025 en comparación con el año 2024; pero si se revisa el ciclo del año cafetero (2024-2025) el panorama es mucho más positivo y se pueden evidenciar periodos de bononza para este sector.

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

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

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

# Verificar longitudes
print(length(fechas_corregidas_var2))
## [1] 156
print(length(tasa_crecimiento_var2))
## [1] 156
print(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 = "#D2B48C", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "#8B4513", size = 0.8, linetype = "dashed") +
  ggtitle("PECAFE: 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)

Se observa un incremento pronunciado del precio externo, explicado por la reducción de la oferta mundial debido al Fenómeno de El Niño, que afectó la producción en Brasil y otros países productores. Colombia se benefició de esta escasez, posicionando su café arábica lavado a precios más altos en los mercados internacionales.

Entre finales de 2014 y Agosto de 2015 El precio externo cae de forma marcada.Esto se explica por la recuperación de la oferta mundial, la fortaleza del dólar y una corrección del mercado internacional tras el repunte previo.Este comportamiento marcó el inicio de un ciclo de ajuste y estabilización que se prolongó hasta 2015.

En 2015 inicia un crecimiento, aunque con tasas negativas, y es hasta el segundo semestre de 2016 que se observan valores positivos en la tasa de crecimiento. Este repunte se mantiene hasta 2017, cuando nuevamente el precio externo pasa a terreno negativo, reflejando la volatilidad del mercado internacional y la recuperación de la oferta en países competidores.

Entre 2020 y 2021 se observa un crecimiento positivo.La pandemia generó disrupciones logísticas y escasez temporal de café en el mercado internacional, lo que elevó los precios externos.

En 2022 los precios se ajustan nuevamente y se inicia un descenso.La recuperación de la producción en Brasil y Vietnam generó un exceso de oferta, reduciendo los precios.Simultáneamente, el Fenómeno de La Niña afectó la calidad del grano en Colombia, limitando su competitividad en mercados premium

Para el año 2023 a 2024 El precio externo alcanza uno de sus niveles más altos de la década.Explicados en la heladas y sequías en Brasil, que redujeron la oferta global tambien hubo mayor demanda de cafés premium y certificados y los costos logísticos elevados, que se trasladaron al precio final.

A inicios de 2025 el precio externo comienza a caer. La oferta mundial se normaliza y los costos logísticos se estabilizan. Además, la depreciación del peso colombiano reduce el incentivo para mantener precios altos en dólares.La tendencia negativa refleja una fase de ajuste del mercado internacional tras los picos de 2023–2024.

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

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

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

# Verificar longitudes
print(length(fechas_corregidas_var3))
## [1] 156
print(length(tasa_crecimiento_var3))
## [1] 156
print(length(tasa_tendencia_var3))
## [1] 156
# 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 = "#D2B48C", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "#8B4513", size = 0.8, linetype = "dashed") +
  ggtitle("XCAF: 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)

Para la tasa de crecimiento de las exportaciones de Cafe, se observa que la variable no se mantien estable,muestra un patrón cíclico con fases de expansión y contracción, típico de un sector dependiente del clima y de la demanda internacional

De 2013 a 2014 Las exportaciones aumentan gracias a la recuperación de la producción nacional tras la crisis de la roya (2008–2011). La renovación de cafetales con la variedad Castillo permitió una mayor oferta exportable. Además, los precios internacionales favorables impulsaron los envíos al exterior y fortalecieron la posición de Colombia en el mercado global.

Entre 2015 y 2016 Se observa una disminución en la tasa de crecimiento, asociada al Fenómeno de El Niño (2015–2016), que afectó la producción y calidad del grano. La menor oferta redujo los volúmenes exportados. También influyó la volatilidad de la TRM, que alteró los márgenes de ganancia de los exportadores y generó incertidumbre en los contratos internacionales.

De 2019 a 2021 se presenta un crecimiento sostenido, posiblemente vinculado a la reposición de inventarios internacionales y al incremento del consumo doméstico de café durante la pandemia. Las restricciones logísticas globales elevaron los precios internacionales y estimularon las exportaciones colombianas, especialmente hacia Estados Unidos y Europa.

De 2022 a 2023 Las exportaciones vuelven a caer debido al Fenómeno de La Niña, que generó exceso de lluvias y afectó la calidad del grano. A esto se sumaron problemas logísticos en puertos y transporte internacional, así como altos costos de insumos y riesgos de seguridad rural, que limitaron la capacidad de cumplir contratos de exportación.

A finales de 2023 y hasta inicios del 2025 Se observa una recuperación estable, impulsada por la demanda de cafés especiales y sostenibles en Europa y Asia. La mejora en la producción nacional y la estabilización de los precios internacionales favorecieron los volúmenes exportados y consolidaron la reputación del café colombiano como producto premium.

A partir de 2025, La tendencia vuelve a caer, explicada por el exceso de lluvias, la menor inversión agronómica y los mayores costos logísticos. La oferta exportable se reduce y los márgenes se estrechan, generando una contracción en la tasa de crecimiento y evidenciando la vulnerabilidad del sector ante factores externos.

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

Modelo ARIMA

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

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

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

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

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

# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 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 oct2025

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

Identificación automática del modelo ARIMA

library(forecast)

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

# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)
## Series: train_ts 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.2548  -0.5399
## s.e.   0.0715   0.0713
## 
## sigma^2 = 32369:  log likelihood = -1083.79
## AIC=2173.59   AICc=2173.74   BIC=2182.89
## 
## Training set error measures:
##                   ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 16.8956 178.2713 141.1356 -0.4770831 13.48926 0.8018914
##                      ACF1
## Training set -0.001962156

Estimación del modelo identificado automatico y validación de Significancia de coeficientes

library(lmtest)

# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.254763   0.071528 -3.5617 0.0003684 ***
## ma2 -0.539876   0.071270 -7.5751  3.59e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajuste del modelo ARIMA(0,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.2548  -0.5399
## s.e.   0.0715   0.0713
## 
## sigma^2 = 32369:  log likelihood = -1083.79
## AIC=2173.59   AICc=2173.74   BIC=2182.89
## 
## Training set error measures:
##                   ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 16.8956 178.2713 141.1356 -0.4770831 13.48926 0.8018914
##                      ACF1
## Training set -0.001962156

Interpretación de Coeficientes

Se puede observar que ambos coeficientes son negativos y altamente significativos (***), lo que significa que el modelo esta corrigiendo sistemáticamente los choques pasados.Adicionalmente se observa que ma2 es más fuerte que ma1 lo cual indica que los errores de dos periodos atrás influyen más en la dinámica actual que los del periodo inmediatamente anterior; esto puede ser coherente o acorde con una serie como la producción de café, donde los efectos de un choque (Clima, Cosechas) no se corrigen de inmediato, si no que se reflejan con más fuerza en los siguientes ciclos.

En cuanto a los resultados de Training Set Error Measures, se puede evidenciar que el MAPE por ejemplo el error porcentual medio es aceptable, de acuerdo con su medición al estar por debajo del 20% puede considerarse útil, si embargo al estar en un 13% aproximadamente el modelo no logra una precisión razonable, aunque no óptima; por otro lado el RMSE es relativamente alto, lo cual indica que aunque el modelo siga la tendencia, no logra capturar bien los picos y caídas abruptas que podrian ser propias o fecuentes en el sector cafetero.

Validación de residuales o errores del modelo

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima_auto)  # Verificar si los residuos son aleatorios y no presentan patrones

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

Con respecto a las gráficas relacionadas, se puede observar que para el primer gráfico de residuales en el tiempo, contiene picos grandes entre -250 y 500, lo que refleja que el modelo no logra capturar del todo los choques abruptos de la serie, esto podria ocurrir en un modelo Arima en el cual se suaviza la serie, pero deja errores grandes cuando hay ciclos o shocks externos.

La gráfica ACF de los residuales, se peude observar que la mayoria de las barras están dentro de las bandas de confianza que representan las lineas azules, sin embargo es observan picos en rezagos como (12, 24 y 36 ) lo que sugiere que podria existir cierta estacionalidad o periocidad en la serie de ARIMA no esta captando.

El histograma muestra una concentración importante en el medio o centrada en cero y menos casos a los extremos pero si un poco largas, lo que refleja que en algunos periodos los errores fueron grandes como choques no capturados;la curva de color naranja se ajusta bien con respecto a su forma de campana, lo cual es algo positivo con respecto a la distribución de los datos.

Pronóstico modelo ARIMA automático dentro de muestra o en el set de prueba

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

# 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 PNCAFE") +
  xlab("Tiempo") + ylab("Miles de sacos de 60 Kg de café verde equivalente")

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

Interpretación modelo automatico (0,1,2): El modelo automático ARIMA(0,1,2) evidencia un comportamiento marcado de sobreestimación y subestimación en la mayoría de los tramos: cuando los datos observados muestran incrementos, el pronóstico proyecta caídas, y viceversa, manteniéndose casi siempre en oposición; la única coincidencia clara se presenta en el punto de quiebre del año 2025, donde ambos reflejan el mismo cambio de dirección. Al final de la gráfica, mientras los valores observados muestran una inclinación descendente en la producción, el pronóstico se mantiene estable, lo que confirma que el modelo es útil para captar la tendencia global, pero insuficiente para representar con precisión la dinámica estacional y los choques propios del sector cafetero, siendo un poco limitante para la toma de decisiones con respecto a la variable de producción.

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

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

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

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

# Mostrar la tabla
print(forecast_table_auto)
##     Tiempo Observado Pronosticado
## 1 2025.750  1208.087     1186.604
## 2 2025.833  1266.070     1152.042
## 3 2025.917  1233.423     1152.042

Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, es decir enero 2026

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

# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_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   1186.604
## 2 2025.833   1152.042
## 3 2025.917   1152.042
## 4 2026.000   1152.042
# 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: 2026 = 1152.04182460378"

Se puede observar que, el modelo proyecta una caida inicial desde 1186 hasta 1152 en el primer tramo, a partir de ahi, los valores se estabilizarn en 1152, sin mostrar variaciones adicionales.Esto refleja que el ARIMA sin estacionalidad tiende a aplanar la serie en el futuro, manteniendo un nivel constante en lugar de anticipar ciclos o fluctuaciones.

Modelo SARIMA automático

Este modelo podria ser una solución o mejora al modelo arima tradicional ya que recoge el efecto estacional de las variables, es recomendable por tanto para datos que si tienen un componente estacional fuerte.

El modelo ajustado en este ejemplo es un SARIMA(0,1,1)(1,0,0)[12], lo que significa:

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

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

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

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

Identificación automática del modelo SARIMA

# Identificación automática modelo SARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)
## Series: train_ts 
## ARIMA(0,1,2)(2,0,0)[12] 
## 
## Coefficients:
##           ma1      ma2    sar1    sar2
##       -0.3560  -0.4240  0.3121  0.3338
## s.e.   0.0776   0.0774  0.0762  0.0787
## 
## sigma^2 = 23617:  log likelihood = -1059.69
## AIC=2129.37   AICc=2129.75   BIC=2144.87

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(0, 1, 1),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(1, 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(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.7226  0.5650
## s.e.   0.1198  0.0715
## 
## sigma^2 = 30375:  log likelihood = -1080.72
## AIC=2167.44   AICc=2167.59   BIC=2176.74
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE    MAPE     MASE      ACF1
## Training set 6.058053 172.6936 134.4539 -1.520063 13.0971 0.763928 0.2570015

Análisis de los coeficientes:

El modelo SARIMA(0,1,1)(1,0,0)[12] refleja un ajuste más sólido al incorporar la estacionalidad anual de la serie, lo que se evidencia en sus resultados numéricos: el error medio (ME ≈ 6.05) es prácticamente nulo, mostrando ausencia de sesgo sistemático; el RMSE (≈ 172.7) y el MAPE (≈ 13.1%) son más bajos que en el ARIMA, confirmando una mayor precisión en las predicciones; y el MPE (≈ -1.52%) indica apenas una leve tendencia a subestimar. Además, el coeficiente estacional autorregresivo (SAR1 = 0.565) demuestra que los valores actuales dependen significativamente de los registros de hace 12 periodos, capturando así los ciclos propios del café.

Validación de residuales del modelo automatico SARIMA

En el correlograma de residuos siguiente se observa que, mejora la correlación de los residuos frente a lso dos modelos anteriores. Sin embargo, al comparar los valores reales VS pronosticados se determina una poca coincidencia. Sigue funcionando mejor el modelo automatico (4,1,2)

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(1,0,0)[12]
## Q* = 65.29, df = 22, p-value = 3.566e-06
## 
## Model df: 2.   Total lags used: 24

Análisis de gráficos:

Para el gráfico uno, se pueden visualizar picos importantes, el escenario no cambia radicalmente frente al modelo ARIMA, sin embargo el modelo de SARIMA se logra absorber mejor la estacionalidad, esto quiere decir que, aunque la magnitud de algunos residuos sigue siendo alta, el modelo evita que se acumulen patrones de error y mantiene los residuos más cercanos a ruido blanco, por lo cual mejora en terminos de estructura, pero no elimina por completo los choques fuertes que son propias del sector cafetero.

Para el gráfico ACF, se puede evidenciar que aunque algunas barras aún sobrepasan las líneas azules de confianza, estas aparecen de manera aislada y sin un patrón repetitivo en los rezagos estacionales; esto significa que no hay una autocorrelación sistemática y que los residuos se comportan en gran medida como ruido blanco.

Para el gráfico de Histograma se puede evidenciar que, la distribución es simétrica y centrada en cero, lo que indica la ausencia de sesgo sistemático. La curva narajan se ajusta bastante bien al histograma, reforzando la idea de que los residuos se comportan de manera aleatoria. Sin embargo, todavía aparecen colas moderadas, lo que significa que existen algunos valores extremos que el modelo no logra capturar del todo.

Pronóstico con el modelo SARIMA dentro del set de prueba-Gráfico líneas

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

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

# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado PNCAFE") +
  xlab("Tiempo") + ylab("Miles de sacos de 60 Kg de café verde equivalente")

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

Análisis del pronóstico:

En este gráfico de Pronóstico vs Observado, se puede evidenciar que, se sigue notando la tendencia de sobrestimacion y subestimación; el pronóstico que corresponde a la linea azul, se ubica por encima de los valores reales en algunos tramos y por debajo en otros. Sin embargo, la diferencia clave es que para este modelo el ajuste resulta más acertado, porque el comportamiento proyectado logra acercarse mejor a la dinámica de la serie, capturando la tendencia general y reflejando con mayor fidelidad los cliclos estacionales. En conclusión aunque persisten discrepancias puntuales, sin embargo el modelo SARIMA consigue reducir la distancia, lo cual arroja un pronóstico más útil para evaluar el comportamiento de la producción de café.

Lo anterior nos permite concluir que el modelo automático de SARIMA (0,1,1) fue el que tuvo mejor desempeño, pese a las especificaciones mencionadas en el párrafo anterior con relación en la sobrestimación y subestimación, este modelo logra capturar la estacionalidad anual presente en la serie del café, reflejando con mayor precisión los ciclos de producción, lo cual hace que el pronóstico sea más coherente y ajustado a la realidad del sector, convirtiendolo en una opción más confiable para el análisis y la planificación.

Pronóstico del modelo automático SARIMA en el set de prueba-Tabla

# 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  1208.087     1342.109
## 2 2025.833  1266.070     1510.135
## 3 2025.917  1233.423     1500.895

Pronóstico del modelo automático SARIMA fuera de muestra, es decir, en enero 2026 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   1342.109
## 2 2025.833   1510.135
## 3 2025.917   1500.895
## 4 2026.000   1275.700
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2026:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2026: 2026 = 1275.6998359887"

Conclusión

Podemos concluir que, el pronóstico del sector cafetero muestra señales de un ligero crecimiento en exportaciones y precios internacionales, lo que refleja la fortaleza de la demanda externa y la reputación del café colombiano en mercados premium. Aunque la producción mantiene un comportamiento volátil y difícil de anticipar con precisión, la extracción de señales de las demás variables permite intuir que los ingresos por exportaciones podrían mejorar en los próximos años, incluso en escenarios de oferta limitada.

Este panorama obliga a reconocer los riesgos: la volatilidad productiva por factores climáticos y plagas, la dependencia de precios internacionales, la exposición cambiaria, la concentración en pocos mercados y el incremento de costos agrícolas que presionan la rentabilidad. Sin embargo, también abre oportunidades estratégicas: aprovechar el repunte de precios internacionales, consolidar la expansión de exportaciones, diferenciarse por calidad en segmentos premium, diversificar hacia Asia y Medio Oriente, y fortalecer la innovación y sostenibilidad con cafés certificados.

Recomendaciones estratégicas del sector

El análisis anterior nos ha permitido concluir que la compañia debería adoptar acciones o una estrategia integral que combine medidas agrícolas, financieras y comerciales, enfocándose en reducir la vulnerabilidad de la producción frente a plagas y fenómenos climáticos mediante la tecnificación del cultivo, la renovación progresiva de cafetales con variedades resistentes y la instalación de sistemas de riego que estabilicen los volúmenes. Al mismo tiempo, conviene implementar coberturas en mercados de futuros y cambiarios para proteger los ingresos ante la volatilidad de precios internacionales, diversificar los destinos de exportación hacia Asia y Medio Oriente para disminuir la dependencia de EE. UU. y Europa, y certificar parte de la producción bajo sellos orgánicos o de comercio justo que permitan acceder a segmentos premium con mayor valor agregado. Estas decisiones concretas fortalecen la resiliencia productiva y financiera, al tiempo que posicionan a la empresa en mercados más rentables y sostenibles.

Indicadores que deberia seguir monitoreando:

1.Producción nacional de café (PNCAFE): volumen total y tasas de crecimiento para anticipar ciclos de caída o recuperación.

2.Precio internacional del café colombiano (PECAFE): evolución en centavos/libra en la Bolsa de Nueva York, clave para ingresos.

3.Exportaciones de café (XCAF): volumen y destinos, identificando concentración o diversificación de mercados.

4.Tasa de cambio COP/USD (TRM): impacto directo en la competitividad y en los ingresos en pesos.

Áreas de la empresa a las que podria servirles esta información:

*Área comercial/ventas: puede usar los pronósticos de precios y exportaciones para definir estrategias de negociación, ajustar márgenes y decidir cuándo cerrar contratos internacionales.

*Área de producción: aprovecha los datos de volatilidad para planificar la cantidad a sembrar, renovar cafetales y organizar recursos según los ciclos de mayor o menor producción.

*Área de comercio exterior: analiza qué mercados son más atractivos en cada periodo, identifica zonas y destinos más provechosos y ajusta la logística de exportación.

*Área Financiera: podría tomar decisiones en relación con la aplicación de coberturas cambiarias o de precios, asegurando estabilidad en los ingresos.