# ── Carga de librerías ────────────────────────────────────────────────────────
library(readxl)        # Lectura de archivos Excel
library(tidyverse)     # Manipulación y visualización de datos
library(lubridate)     # Manejo de fechas
library(forecast)      # Modelos ARIMA / SARIMA
library(tseries)       # Pruebas de raíz unitaria (ADF)
library(zoo)           # Objetos de series de tiempo irregulares
library(ggplot2)       # Gráficos
library(gridExtra)     # Composición de gráficos
library(scales)        # Formato de ejes
library(knitr)         # Tablas en HTML
library(kableExtra)    # Tablas enriquecidas
library(seasonal)      # Descomposición STL / X-13
library(FinTS)         # Prueba ARCH
library(Metrics)       # Métricas de error (RMSE, MAE)

1 Contexto: Sector, Empresa y Variables

1.1 El Sector Agroindustrial Azucarero del Valle del Cauca

El sector agroindustrial del azúcar constituye uno de los pilares históricos de la economía del Valle del Cauca y del suroccidente colombiano. Colombia es el único país del mundo que cosecha caña de azúcar durante los doce meses del año, gracias a las condiciones geoclimáticas únicas del valle geográfico del río Cauca, ubicado entre 900 y 1.100 metros sobre el nivel del mar, con temperatura promedio de 24 °C, alta luminosidad solar y disponibilidad hídrica de los ríos que descienden de las cordilleras Occidental y Central.

La cadena productiva se extiende desde el cultivo de la caña (Saccharum officinarum), pasando por la molienda en los ingenios, hasta la producción de azúcar crudo, azúcar refinada, mieles, alcohol carburante y cogeneración de energía eléctrica. Existen actualmente trece ingenios en la región, siendo los más representativos Ingenio del Cauca (Incauca), Manuelita, Providencia, Mayagüez y Risaralda, entre otros.

1.2 La Empresa: Ingenio Providencia S.A. (Empresa de Referencia)

Para efectos de este análisis, se toma como empresa de referencia a Ingenio Providencia S.A., uno de los conglomerados azucareros más grandes de Colombia, con presencia en producción de azúcar, alcohol carburante, energía y derivados alimenticios. La empresa es representativa del sector dado que sus resultados reflejan fielmente las dinámicas macroeconómicas, productivas y de comercio exterior que afectan al conjunto de la industria azucarera vallecaucana.

El ingenio procesa anualmente más de dos millones de toneladas de caña, produce azúcar blanca y cruda, y participa activamente en el mercado de exportaciones agroindustriales del departamento. Su gestión estratégica debe responder tanto a ciclos productivos estacionales propios del cultivo, como a la volatilidad de los mercados internacionales y las condiciones macroeconómicas domésticas.

1.3 Justificación y Descripción de las Tres Variables

Las tres variables seleccionadas capturan dimensiones fundamentales e interconectadas del ciclo agroindustrial: la capacidad productiva, la transformación industrial y la inserción en los mercados externos.

vars_df <- data.frame(
  `N°`       = 1:3,
  Variable   = c("Molienda de Caña de Azúcar",
                 "Exportaciones Totales del Valle del Cauca",
                 "Producción de Azúcar"),
  Acrónimo   = c("CAN", "X_V", "AZUCAR"),
  Unidad     = c("Toneladas", "Dólares CIF", "Toneladas"),
  Geografía  = c("Colombia", "Valle del Cauca", "Colombia"),
  Fuente     = c("DANE", "DANE", "ASOCAÑA / DANE"),
  Justificación = c(
    "Es el insumo primario y el primer eslabón de la cadena. Su nivel refleja la disponibilidad de materia prima, las decisiones de cosecha, el estado fitosanitario del cultivo y las condiciones climáticas. Una mayor molienda es condición necesaria —aunque no suficiente— para una mayor producción de azúcar.",
    "Captura la competitividad del departamento en el mercado mundial. El Valle del Cauca exporta principalmente azúcares, preparaciones alimenticias, productos químicos y papel. La variable refleja el acceso a divisas, el impacto de la tasa de cambio y la demanda internacional de los productos agroindustriales.",
    "Es el output final del proceso industrial y el principal producto generador de ingresos para los ingenios. Su nivel determina directamente los ingresos operativos y la rentabilidad empresarial. Su comportamiento respecto a la molienda revela la eficiencia extractiva (rendimiento caña-azúcar)."
  ),
  check.names = FALSE
)

kable(vars_df[, 1:6], caption = "Tabla 1. Catálogo de Variables Seleccionadas",
      align = c("c","l","c","l","c","l")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE, font_size = 12) %>%
  column_spec(1, bold = TRUE, width = "3em") %>%
  column_spec(2, bold = TRUE, color = "#1a3c5e") %>%
  column_spec(3, monospace = TRUE, color = "#e67e22")
Tabla 1. Catálogo de Variables Seleccionadas
N° Variable Acrónimo Unidad Geografía Fuente
1 Molienda de Caña de Azúcar CAN Toneladas Colombia DANE
2 Exportaciones Totales del Valle del Cauca X_V Dólares CIF Valle del Cauca DANE
3 Producción de Azúcar AZUCAR Toneladas Colombia ASOCAÑA / DANE

2 Preparación y Descripción de los Datos

2.1 Carga y Estructuración

# ── Lectura del archivo Excel ─────────────────────────────────────────────────

df_raw <- read_excel("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"))

# ── Selección de variables de interés ────────────────────────────────────────

df <- df_raw %>%
  select(FECHA, CAN, AZUCAR, X_V) %>%
  mutate(FECHA = as.Date(FECHA)) %>%
  arrange(FECHA)

# ── Vista preliminar ──────────────────────────────────────────────────────────
kable(head(df, 12),
      caption = "Tabla 2. Primeras 12 observaciones — Sector Agroindustrial del Azúcar en el Valle del Cauca",
      col.names = c("Fecha", "Molienda de Caña de azúcar (Toneladas)", "Producción de azúcar (Toneladas)", "Exportaciones totales del Valle (Dólares CIF)"),
      format.args = list(big.mark = ",", digits = 2),
      align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 12)
Tabla 2. Primeras 12 observaciones — Sector Agroindustrial del Azúcar en el Valle del Cauca
Fecha Molienda de Caña de azúcar (Toneladas) Producción de azúcar (Toneladas) Exportaciones totales del Valle (Dólares CIF)
2012-01-01 1,685,584 148,483 1.6e+08
2012-02-01 1,973,654 192,196 1.8e+08
2012-03-01 2,083,470 202,408 1.9e+08
2012-04-01 1,406,868 134,532 1.6e+08
2012-05-01 1,233,631 107,659 2.0e+08
2012-06-01 1,997,318 206,134 1.7e+08
2012-07-01 2,107,652 222,915 1.9e+08
2012-08-01 2,075,408 226,868 2.0e+08
2012-09-01 2,005,429 225,330 1.9e+08
2012-10-01 1,700,440 169,420 2.2e+08
2012-11-01 1,144,776 103,433 1.8e+08
2012-12-01 1,409,399 138,275 1.9e+08

2.2 Estadísticas Descriptivas

# ── Tabla de estadísticas descriptivas ───────────────────────────────────────
desc_stats <- df %>%
  select(-FECHA) %>%
  summarise(across(everything(), list(
    N         = ~n(),
    Media     = ~mean(., na.rm = TRUE),
    Mediana   = ~median(., na.rm = TRUE),
    Des.Est.  = ~sd(., na.rm = TRUE),
    Mínimo    = ~min(., na.rm = TRUE),
    Máximo    = ~max(., na.rm = TRUE),
    CV        = ~(sd(., na.rm = TRUE)/mean(., na.rm = TRUE))*100
  ))) %>%
  pivot_longer(everything(),
               names_to = c("Variable", "Estadístico"),
               names_sep = "_(?=[^_]+$)") %>%
  pivot_wider(names_from = "Estadístico", values_from = "value")

kable(desc_stats,
      caption = "Tabla 3. Estadísticas Descriptivas (2012–2025)",
      digits = 2,
      format.args = list(big.mark = ","),
      align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE, font_size = 12)
Tabla 3. Estadísticas Descriptivas (2012–2025)
Variable N Media Mediana Des.Est. Mínimo Máximo CV
CAN 168 1,920,333 1,997,577.7 396,510.60 138,148.00 2,535,470.3 20.65
AZUCAR 168 179,856 183,641.6 42,611.88 11,041.41 263,433.5 23.69
X_V 168 157,545,398 157,166,984.6 24,349,256.05 49,767,731.10 222,055,829.3 15.46

2.3 Construcción de Series de Tiempo

# ── Series de tiempo mensuales (frecuencia = 12) ──────────────────────────────
ts_CAN <- ts(df$CAN, start = c(2012, 1), frequency = 12)
ts_AZUCAR  <- ts(df$AZUCAR,  start = c(2012, 1), frequency = 12)
ts_X_V    <- ts(df$X_V,    start = c(2012, 1), frequency = 12)

3 Extracción de Señales

3.1 Evolución Histórica de las Series (2012–2025)

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)

# Definir paleta de colores
PAL <- c(
  "Molienda de Caña (ton)" = "#1b9e77",
  "Exportaciones Valle (USD CIF)" = "#d95f02",
  "Producción de Azúcar (ton)" = "#7570b3"
)

# Tema (por si no lo tienes)
theme_sector <- function() {
  theme_minimal()
}

# Convertir a formato largo
datos_long <- df %>%
  pivot_longer(cols = c(CAN, X_V, AZUCAR),
               names_to = "Variable", values_to = "Valor") %>%
  mutate(Variable = recode(Variable,
    CAN    = "Molienda de Caña (ton)",
    X_V    = "Exportaciones Valle (USD CIF)",
    AZUCAR = "Producción de Azúcar (ton)"))

# Gráfico
ggplot(datos_long, aes(x = FECHA, y = Valor, colour = Variable)) +
  geom_line(linewidth = 0.8, alpha = 0.9) +
  geom_smooth(method = "loess", se = FALSE, linewidth = 0.4,
              linetype = "dashed", alpha = 0.5) +
  facet_wrap(~Variable, scales = "free_y", ncol = 1) +
  scale_colour_manual(values = PAL) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(labels = label_comma()) +
  labs(
    title    = "Figura 1. Evolución Histórica de las Variables del Sector Azucarero",
    subtitle = "Series mensuales 2012–2025 | Línea discontinua: tendencia LOESS",
    x = NULL, y = NULL,
    caption  = "Fuente: DANE – ASOCAÑA | Elaboración propia"
  ) +
  theme_sector() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))

4 Análisis Exploratorio: Evolución Histórica

4.1 Molienda de Caña de azúcar (CAN)

df_CAN <- as.data.frame(window(ts_CAN, start = c(2012, 1)))
df_CAN$Index <- time(window(ts_CAN, start = c(2012, 1)))

autoplot(window(ts_CAN, start = c(2012, 1)), colour = "#4E6D4E") +
  geom_smooth(data = df_CAN,
              aes(x = Index, y = df_CAN[,1]),
              method   = "loess",
              colour   = "#2C3E2D",
              linetype = "dashed",
              se       = FALSE,
              linewidth = 0.8) +
  labs(
    title    = "Molienda de Caña de azúcar",
    subtitle = "Toneladas — Ene 2012 a Dic 2025",
    x        = "Período",
    y        = "Toneladas",
    caption  = "Fuente: DANE"
  ) +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))
Figura 1. Molienda de Caña de azúcar (2012–2025)

Figura 1. Molienda de Caña de azúcar (2012–2025)

La producción de café en Colombia exhibe una dinámica marcadamente cíclica con patrones estacionales pronunciados asociados a las dos cosechas anuales. Durante el período analizado se observa una tendencia positiva moderada desde 2012 hasta aproximadamente 2017–2018, período en el que la producción alcanzó niveles cercanos a los 1.700–1.800 miles de sacos mensuales, impulsada por los programas de renovación de cafetales de la Federación Nacional de Cafeteros. Posteriormente, se registra una fase de corrección y mayor volatilidad, con impactos derivados del fenómeno de La Niña (exceso de lluvias) y las perturbaciones de la pandemia de COVID-19 en 2020–2021. Hacia el período 2022–2025, la producción muestra señales de recuperación, aunque con oscilaciones más amplias que en el período inicial.

4.2 Producción de azúcar (AZUCAR)

library(ggplot2)
library(scales)

# Convertir serie a dataframe correctamente
ts_data <- ts_AZUCAR / 1e6

df_azucar <- data.frame(
  Index = time(ts_data),
  Valor = as.numeric(ts_data)
)

# Gráfico
autoplot(ts_data, colour = "#8B4513") +
  geom_smooth(data = df_azucar,
              aes(x = Index, y = Valor),
              method    = "loess",
              colour    = "#5C2D0A",
              linetype  = "dashed",
              se        = FALSE,
              linewidth = 0.8) +
  labs(
    title    = "Producción de azúcar",
    subtitle = "Toneladas — Ene 2012 a Dic 2025",
    x        = "Período",
    y        = "Millones de toneladas",
    caption  = "Fuente: DANE"
  ) +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))
Figura 2. Producción de azúcar (2012–2025)

Figura 2. Producción de azúcar (2012–2025)

El precio interno del café colombiano presenta una tendencia alcista sostenida y pronunciada a lo largo del período analizado. Partiendo de niveles cercanos a los 400–900 mil pesos por carga en 2012–2015, el precio experimenta incrementos significativos hacia 2021–2022, alcanzando máximos históricos superiores a los 3 millones de pesos por carga en el período 2024–2025. Este comportamiento refleja la confluencia de múltiples factores: la depreciación acumulada del peso colombiano frente al dólar (la TRM pasó de ~1.800 a más de 4.000 pesos por dólar), el aumento del precio internacional del café en la bolsa de Nueva York (afectado por problemas de oferta en Brasil y Vietnam), y la inflación doméstica que incrementa los costos de producción internos. El precio interno tiene un componente estacional moderado, pero su volatilidad está dominada por el componente irregular y la tendencia de largo plazo.

4.3 Exportaciones totales del Valle (X_V)

autoplot(ts_X_V / 1e3, colour = "#D4A017") +
  geom_smooth(method = "loess",
              colour = "#8B6914",
              linetype = "dashed",
              se = FALSE,
              linewidth = 0.8) +
  labs(
    title    = "Exportaciones totales del Valle",
    subtitle = "Dólares CIF — Ene 2012 a Dic 2025",
    x        = "Período",
    y        = "Dólares CIF (miles)",
    caption  = "Fuente: DANE — Federación Nacional de Cafeteros de Colombia"
  ) +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))
Figura 3. Exportaciones mensuales de café colombiano (2012–2025)

Figura 3. Exportaciones mensuales de café colombiano (2012–2025)

Las exportaciones de café muestran un patrón de crecimiento notable a lo largo de los 14 años analizados, pasando de un promedio mensual de aproximadamente 115–200 millones de dólares en 2012 a niveles que superan los 400–590 millones de dólares en el período 2024–2025. Esta trayectoria ascendente responde tanto al aumento del precio internacional como al posicionamiento diferenciado del café colombiano en mercados de especialidad. Se identifican fuertes oscilaciones mensuales vinculadas al ciclo productivo y a los patrones de demanda de los mercados de destino (mayor demanda en los meses fríos del hemisferio norte). El año 2020 muestra una contracción temporal asociada al cierre de puertos y las disrupciones logísticas de la pandemia de COVID-19.

5 Extracción de Señales

En esta sección se aplica la descomposición STL (Seasonal and Trend decomposition using Loess), método robusto que permite extraer de manera independiente los componentes de tendencia, estacionalidad y componente irregular (residuo) de cada serie de tiempo.

5.1 Descomposición STL — MOLIENDA DE CAÑA DE AZÚCAR (CAN)

library(plotly)

stl_AZUCAR <- stl(ts_AZUCAR, s.window = "periodic", robust = TRUE)

df_stl_azucar <- data.frame(
  fecha      = as.Date(time(ts_AZUCAR)),
  original   = as.numeric(ts_AZUCAR),
  estacional = as.numeric(stl_AZUCAR$time.series[, "seasonal"]),
  tendencia  = as.numeric(stl_AZUCAR$time.series[, "trend"]),
  residuo    = as.numeric(stl_AZUCAR$time.series[, "remainder"])
)

colores <- c(
  original   = "#8B4513",
  tendencia  = "#5C2D0A",
  estacional = "#D2956A",
  residuo    = "#C0392B"
)

# ── Paneles SIN layout individual ─────────────────────────
p1 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~original,
              type = "scatter", mode = "lines",
              name = "Original",
              line = list(color = colores["original"], width = 1.5))

p2 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~estacional,
              type = "scatter", mode = "lines",
              name = "Estacionalidad",
              line = list(color = colores["estacional"], width = 1.5))

p3 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~tendencia,
              type = "scatter", mode = "lines",
              name = "Tendencia",
              line = list(color = colores["tendencia"], width = 2))

p4 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~residuo,
              type = "bar", name = "Residuo",
              marker = list(
                color   = ifelse(df_stl_azucar$residuo >= 0,
                                 colores["tendencia"],
                                 colores["residuo"]),
                opacity = 0.7
              ))

# ── Subplot con títulos de ejes en layout global ───────────
subplot(p1, p2, p3, p4,
        nrows = 4, shareX = TRUE, titleY = TRUE,
        heights = c(0.30, 0.25, 0.25, 0.20)) %>%
  layout(
    title = list(
      text = "<b>Descomposición STL: Producción de Azúcar</b><br>
              <sup>Fuente: Federación Nacional de Cafeteros / Elaboración propia</sup>",
      font = list(size = 15)
    ),
    yaxis  = list(title = "Original"),
    yaxis2 = list(title = "Estacional"),
    yaxis3 = list(title = "Tendencia"),
    yaxis4 = list(title = "Residuo"),
    hovermode     = "x unified",
    paper_bgcolor = "#FAFAFA",
    plot_bgcolor  = "#FAFAFA",
    legend = list(orientation = "h", x = 0.3, y = -0.05)
  )

Figura 4. Descomposición STL — Producción Nacional de AZUCAR (CAN)

5.1.1 Interpretación — CAN

Tendencia: El componente tendencial de la producción de café muestra una dinámica en tres fases claramente diferenciadas. La primera fase (2012–2019) refleja una tendencia creciente moderada, que pasó de aproximadamente 850 miles de sacos mensuales a niveles superiores a los 1.100–1.200 miles. Esta expansión fue impulsada por los programas de renovación de cafetales ejecutados por la Federación Nacional de Cafeteros desde 2008 y que alcanzaron su madurez productiva entre 2015 y 2019. La segunda fase (2019–2021) evidencia una corrección a la baja atribuible a la combinación del fenómeno climático de La Niña, que generó exceso de precipitaciones en las zonas cafeteras, y las restricciones logísticas derivadas de la pandemia de COVID-19. La tercera fase (2022–2025) muestra una tenue recuperación hacia los 1.100–1.150 miles de sacos, aunque sin alcanzar los picos históricos observados en el período previo.

Estacionalidad: La producción de café en Colombia exhibe un patrón estacional bien definido y estadísticamente significativo, lo que confirma la existencia de dos cosechas anuales: la cosecha principal (octubre–febrero), que concentra aproximadamente el 60–65% de la producción anual, y la cosecha mitaca (abril–julio), de menor magnitud. Los meses de diciembre y enero registran consistentemente los valores más altos del componente estacional, mientras que agosto y septiembre corresponden al período de transición de menor producción entre las dos cosechas. Este patrón es biológicamente determinado por el ciclo de floración y desarrollo del grano del cafeto, siendo relativamente estable en el tiempo.

Componente Irregular: El residuo presenta una varianza moderada a lo largo del período, con picos de alta volatilidad en 2020 (pandemia) y en los años de fenómenos climáticos extremos (2015–2016 con El Niño). Estos shocks irregulares son difícilmente predecibles mediante modelos estadísticos y representan el principal factor de incertidumbre para la planificación de la empresa.

5.2 Descomposición STL — Precio Interno del AZUCAR (AZUCAR)

library(plotly)
library(dplyr)

stl_AZUCAR <- stl(ts_AZUCAR, s.window = "periodic", robust = TRUE)

# Extraer componentes
df_stl_azucar <- data.frame(
  fecha      = as.Date(time(ts_AZUCAR)),
  original   = as.numeric(ts_AZUCAR),
  estacional = as.numeric(stl_AZUCAR$time.series[, "seasonal"]),
  tendencia  = as.numeric(stl_AZUCAR$time.series[, "trend"]),
  residuo    = as.numeric(stl_AZUCAR$time.series[, "remainder"])
)

# Paleta de colores
colores <- c(
  original   = "#8B4513",
  tendencia  = "#5C2D0A",
  estacional = "#D2956A",
  residuo    = "#C0392B"
)

# ── Paneles ────────────────────────────────────────────────
p1 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~original, type = "scatter",
              mode = "lines", name = "Original",
              line = list(color = colores["original"], width = 1.5)) %>%
  layout(yaxis = list(title = "Original"))

p2 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~estacional, type = "scatter",
              mode = "lines", name = "Estacionalidad",
              line = list(color = colores["estacional"], width = 1.5)) %>%
  layout(yaxis = list(title = "Estacional"))

p3 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~tendencia, type = "scatter",
              mode = "lines", name = "Tendencia",
              line = list(color = colores["tendencia"], width = 2)) %>%
  layout(yaxis = list(title = "Tendencia"))

p4 <- plot_ly(df_stl_azucar, x = ~fecha, y = ~residuo, type = "bar",
              name = "Residuo",
              marker = list(color = ifelse(df_stl_azucar$residuo >= 0,
                                           colores["tendencia"],
                                           colores["residuo"]),
                            opacity = 0.7)) %>%
  layout(yaxis = list(title = "Residuo"))

# ── Unir paneles ───────────────────────────────────────────
subplot(p1, p2, p3, p4,
        nrows = 4, shareX = TRUE, titleY = TRUE,
        heights = c(0.30, 0.25, 0.25, 0.20)) %>%
  layout(
    title = list(
      text = "<b>Descomposición STL: Producción de Azúcar</b><br>
              <sup>Fuente: DANE / Elaboración propia</sup>",
      font = list(size = 15)
    ),
    hovermode     = "x unified",
    paper_bgcolor = "#FAFAFA",
    plot_bgcolor  = "#FAFAFA",
    legend = list(orientation = "h", x = 0.3, y = -0.05)
  )

Figura 5. Descomposición STL — Precio Interno Base de Compra del Café (AZUCAR)

5.2.1 Interpretación — AZUCAR

Tendencia: El componente de tendencia del precio interno del café es el más dramático de las tres variables analizadas. La tendencia exhibe un crecimiento prácticamente continuo desde el inicio del período, con una fase de aceleración marcada a partir de 2020–2021. Durante 2012–2019, la tendencia fue relativamente plana con oscilaciones moderadas, reflejando un período de precios internacionales relativamente contenidos y una tasa de cambio que, si bien estaba depreciándose, lo hacía de forma gradual. A partir de 2020, la aceleración del proceso de depreciación del peso colombiano (impulsada por la caída del precio del petróleo, la pandemia, y los choques fiscales del país) combinada con el aumento del precio internacional del café generó una espiral alcista que llevó el precio interno a niveles de 2,5–3,1 millones de pesos por carga al cierre de 2025.

Estacionalidad: A diferencia de la producción, el precio interno del café exhibe un componente estacional de menor magnitud relativa, aunque estadísticamente presente. El patrón estacional refleja la dinámica de oferta y demanda: en los meses de cosecha alta (noviembre–enero), la mayor oferta de café pergamino tiende a moderar los precios de compra; mientras que en los meses de menor producción (agosto–septiembre), la escasez relativa eleva el precio. Sin embargo, este efecto estacional es ampliamente dominado por la volatilidad de la tendencia y el componente irregular, lo que significa que el ciclo productivo tiene un efecto secundario sobre el precio doméstico.

Componente Irregular: El residuo del precio interno presenta oscilaciones de alta frecuencia e intensidad, especialmente a partir de 2021–2022. Esto refleja la sensibilidad de esta variable a noticias externas (variaciones en la bolsa de Nueva York, clima en Brasil, tensiones geopolíticas) que se transmiten al precio interno con una frecuencia prácticamente diaria, pero que en la frecuencia mensual se manifiestan como shocks irregulares de amplitud creciente.

5.3 Descomposición STL — Exportaciones totales del Valle (X_V)

library(plotly)
library(dplyr)

stl_X_V <- stl(ts_X_V, s.window = "periodic", robust = TRUE)

# Extraer componentes
df_stl_xv <- data.frame(
  fecha      = as.Date(time(ts_X_V)),
  original   = as.numeric(ts_X_V),
  estacional = as.numeric(stl_X_V$time.series[, "seasonal"]),
  tendencia  = as.numeric(stl_X_V$time.series[, "trend"]),
  residuo    = as.numeric(stl_X_V$time.series[, "remainder"])
)

# Paleta de colores — dorados/Exportaciones totales del Valle
colores <- c(
  original   = "#D4A017",
  tendencia  = "#8B6914",
  estacional = "#F0C84A",
  residuo    = "#C0392B"
)

# ── Paneles ────────────────────────────────────────────────
p1 <- plot_ly(df_stl_xv, x = ~fecha, y = ~original, type = "scatter",
              mode = "lines", name = "Original",
              line = list(color = colores["original"], width = 1.5)) %>%
  layout(yaxis = list(title = "Original"))

p2 <- plot_ly(df_stl_xv, x = ~fecha, y = ~estacional, type = "scatter",
              mode = "lines", name = "Estacionalidad",
              line = list(color = colores["estacional"], width = 1.5)) %>%
  layout(yaxis = list(title = "Estacional"))

p3 <- plot_ly(df_stl_xv, x = ~fecha, y = ~tendencia, type = "scatter",
              mode = "lines", name = "Tendencia",
              line = list(color = colores["tendencia"], width = 2)) %>%
  layout(yaxis = list(title = "Tendencia"))

p4 <- plot_ly(df_stl_xv, x = ~fecha, y = ~residuo, type = "bar",
              name = "Residuo",
              marker = list(color = ifelse(df_stl_xv$residuo >= 0,
                                           colores["tendencia"],
                                           colores["residuo"]),
                            opacity = 0.7)) %>%
  layout(yaxis = list(title = "Residuo"))

# ── Unir paneles ───────────────────────────────────────────
subplot(p1, p2, p3, p4,
        nrows = 4, shareX = TRUE, titleY = TRUE,
        heights = c(0.30, 0.25, 0.25, 0.20)) %>%
  layout(
    title = list(
      text = "<b>Descomposición STL: Exportaciones totales del Valle </b><br>
              <sup>Fuente: DANE — / Elaboración propia</sup>",
      font = list(size = 15)
    ),
    hovermode     = "x unified",
    paper_bgcolor = "#FAFAFA",
    plot_bgcolor  = "#FAFAFA",
    legend = list(orientation = "h", x = 0.3, y = -0.05)
  )

Figura 6. Descomposición STL — Exportaciones totales del Valle (X_V)

5.3.1 Interpretación — X_V

Tendencia: Las exportaciones de café exhiben una tendencia ascendente robusta y sostenida a lo largo de todo el período analizado. Esta tendencia es el resultado de la interacción entre tres fuerzas: el incremento del precio internacional del café (que eleva el valor en dólares por unidad exportada), la recuperación gradual de la producción nacional (que amplía el volumen disponible para exportar), y la depreciación del peso colombiano (que hace más rentable exportar en relación con vender en el mercado doméstico). Los niveles de exportación pasaron de un promedio mensual de aproximadamente 150–200 millones de dólares en 2012 a un promedio superior a los 400 millones de dólares en el período 2023–2025, lo que representa un crecimiento acumulado del orden del 150–200% en términos nominales de dólares.

Estacionalidad: Las exportaciones presentan un patrón estacional moderado pero visible, que se alinea parcialmente con el ciclo de cosecha. Los meses de diciembre, enero y marzo tienden a registrar picos exportadores, en correspondencia con la disponibilidad de café proveniente de la cosecha principal. Los meses de agosto y septiembre, en contraste, muestran valores más bajos del componente estacional, coincidiendo con el período de menor disponibilidad de grano. Este patrón es importante para CafExport desde el punto de vista de la logística de exportación y la negociación de contratos con clientes internacionales.

Componente Irregular: El componente residual de las exportaciones es marcadamente más volátil que el de la producción, lo que refleja la influencia de factores externos (disrupciones logísticas portuarias, variaciones en la demanda de los mercados de destino, negociaciones puntuales de grandes contratos) que afectan los valores de exportación en períodos específicos. El período 2020 muestra el mayor shock negativo del período (pandemia), seguido de una rápida recuperación en 2021.

6 Tasa de Crecimiento Anual

library(reactable)
library(htmltools)
library(dplyr)

df_anual <- df %>%
  mutate(anio = as.integer(format(as.Date(FECHA), "%Y"))) %>%
  group_by(anio) %>%
  summarise(
    CAN_anual    = sum(CAN,     na.rm = TRUE),
    AZUCAR_anual = mean(AZUCAR, na.rm = TRUE),
    X_V_anual    = sum(X_V,     na.rm = TRUE)
  ) %>%
  mutate(
    TCA_CAN    = (CAN_anual    / lag(CAN_anual)    - 1) * 100,
    TCA_AZUCAR = (AZUCAR_anual / lag(AZUCAR_anual) - 1) * 100,
    TCA_X_V    = (X_V_anual    / lag(X_V_anual)    - 1) * 100
  )

# ── Función para colorear TCA ──────────────────────────────
tca_cell <- function(value) {
  if (is.na(value)) return(span("—", style = "color: gray;"))
  color  <- if (value > 0) "#1a7a1a" else "#cc0000"
  flecha <- if (value > 0) "▲" else "▼"
  span(
    style = paste0("color:", color, "; font-weight:bold;"),
    paste0(flecha, " ", formatC(value, format = "f", digits = 2), "%")
  )
}

# ── Función para barras de fondo ───────────────────────────
bar_cell <- function(value, max_val, color) {
  pct <- round(value / max_val * 100)
  div(
    style = paste0(
      "background: linear-gradient(90deg, ", color, "33 ", pct,
      "%, transparent ", pct, "%);",
      "padding: 4px 8px; border-radius: 4px;"
    ),
    formatC(value, format = "f", digits = 0, big.mark = ",")
  )
}

max_CAN    <- max(df_anual$CAN_anual,    na.rm = TRUE)
max_AZUCAR <- max(df_anual$AZUCAR_anual, na.rm = TRUE)
max_XV     <- max(df_anual$X_V_anual,    na.rm = TRUE)

# ── Tabla reactable ────────────────────────────────────────
tabla <- reactable(
  df_anual,
  striped         = TRUE,
  highlight       = TRUE,
  bordered        = TRUE,
  compact         = FALSE,
  defaultPageSize = 20,
  style           = list(fontFamily = "Arial, sans-serif", fontSize = "13px"),

  columns = list(
    anio = colDef(
      name  = "Año",
      align = "center",
      width = 60,
      style = list(fontWeight = "bold", color = "#333")
    ),
    CAN_anual = colDef(
      name  = "Molienda de Caña (ton)",
      align = "center",
      cell  = function(value) bar_cell(value, max_CAN, "#4E6D4E")
    ),
    AZUCAR_anual = colDef(
      name  = "Producción Azúcar (ton)",
      align = "center",
      cell  = function(value) bar_cell(value, max_AZUCAR, "#8B4513")
    ),
    X_V_anual = colDef(
      name  = "Exportaciones Valle (USD)",
      align = "center",
      cell  = function(value) bar_cell(value, max_XV, "#D4A017")
    ),
    TCA_CAN = colDef(
      name  = "TCA Caña (%)",
      align = "center",
      cell  = function(value) tca_cell(value)
    ),
    TCA_AZUCAR = colDef(
      name  = "TCA Azúcar (%)",
      align = "center",
      cell  = function(value) tca_cell(value)
    ),
    TCA_X_V = colDef(
      name  = "TCA Exportaciones (%)",
      align = "center",
      cell  = function(value) tca_cell(value)
    )
  )
)

# ── Título y tabla final ───────────────────────────────────
htmltools::browsable(
  tagList(
    tags$h3(
      style = "font-family:Arial; text-align:center; color:#2C3E2D; margin-bottom:4px;",
      "Tabla 3. Tasa de Crecimiento Anual (%) — Caña, Azúcar y Exportaciones"
    ),
    tags$p(
      style = "font-family:Arial; text-align:center; font-size:12px; color:gray;",
      "Fuente: DANE — Federación Nacional de Cafeteros / Elaboración propia"
    ),
    tabla
  )
)

Tabla 3. Tasa de Crecimiento Anual (%) — Caña, Azúcar y Exportaciones

Fuente: DANE — Federación Nacional de Cafeteros / Elaboración propia

df_anual_long <- df_anual %>%
  filter(!is.na(TCA_CAN)) %>%
  select(anio, TCA_CAN, TCA_AZUCAR, TCA_X_V) %>%
  pivot_longer(-anio, names_to = "Variable", values_to = "TCA") %>%
  mutate(Variable = recode(Variable,
    TCA_CAN = "Producción (CAN)",
    TCA_AZUCAR = "Precio Interno (AZUCAR)",
    TCA_X_V   = "Exportaciones (X_V)"
  ))

ggplot(df_anual_long, aes(x = factor(anio), y = TCA, fill = TCA > 0)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~Variable, scales = "free_y", ncol = 1) +
  geom_hline(yintercept = 0, colour = "black", linewidth = 0.5) +
  scale_fill_manual(values = c("TRUE" = "#27AE60", "FALSE" = "#C0392B")) +
  labs(
    title    = "Tasa de Crecimiento Anual por Variable Cafetera",
    subtitle = "Verde = crecimiento positivo | Rojo = contracción",
    x        = "Año",
    y        = "TCA (%)",
    caption  = "Fuente: Federación Nacional de Cafeteros / DANE — Elaboración propia"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    strip.text   = element_text(face = "bold", size = 11),
    plot.title   = element_text(face = "bold"),
    axis.text.x  = element_text(angle = 45, hjust = 1)
  )
Figura 7. Tasa de Crecimiento Anual (%) — Variables Cafeteras

Figura 7. Tasa de Crecimiento Anual (%) — Variables Cafeteras

7 Relación entre Variables: Señales Conjuntas

7.1 Correlación y Señales del Sector

cor_matrix <- df %>%
  select(CAN, AZUCAR, X_V) %>%
  cor(use = "complete.obs")

kable(round(cor_matrix, 4),
      caption = "Tabla 4. Matriz de Correlación de Pearson — Variables Cafeteras",
      align   = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 13)
Tabla 4. Matriz de Correlación de Pearson — Variables Cafeteras
CAN AZUCAR X_V
CAN 1.0000 0.9606 0.0862
AZUCAR 0.9606 1.0000 0.1810
X_V 0.0862 0.1810 1.0000 
p1 <- ggplot(df, aes(x = CAN, y = X_V/1000)) +
  geom_point(alpha = 0.5, colour = "#27AE60") +
  geom_smooth(method = "lm", colour = "#1A5C30") +
  labs(title = "Producción vs Exportaciones", x = "Miles de sacos", y = "Millones USD") +
  theme_minimal()

p2 <- ggplot(df, aes(x = AZUCAR/1e6, y = X_V/1000)) +
  geom_point(alpha = 0.5, colour = "#8B4513") +
  geom_smooth(method = "lm", colour = "#4A1E05") +
  labs(title = "Precio vs Exportaciones", x = "Millones COP/carga", y = "Millones USD") +
  theme_minimal()

p3 <- ggplot(df, aes(x = CAN, y = AZUCAR/1e6)) +
  geom_point(alpha = 0.5, colour = "#D4A017") +
  geom_smooth(method = "lm", colour = "#7A5700") +
  labs(title = "Producción vs Precio", x = "Miles de sacos", y = "Millones COP/carga") +
  theme_minimal()

grid.arrange(p1, p2, p3, ncol = 2)
Figura 8. Relación entre variables cafeteras — Diagramas de dispersión

Figura 8. Relación entre variables cafeteras — Diagramas de dispersión 

# ── Extracción de componentes ajustados por estacionalidad ─

# CAÑA
st_CAN <- stl(ts_CAN, s.window = "periodic")

# AZÚCAR
stl_AZUCAR <- stl(ts_AZUCAR, s.window = "periodic")

# EXPORTACIONES
stl_X_V <- stl(ts_X_V, s.window = "periodic")


# ── Extracción de series ajustadas ─────────────────────────

ts_CAN_ajustada <- ts_CAN - st_CAN$time.series[, "seasonal"]
ts_AZUCAR_ajustada <- ts_AZUCAR - stl_AZUCAR$time.series[, "seasonal"]
ts_X_V_ajustada <- ts_X_V - stl_X_V$time.series[, "seasonal"]


# ── Verificación ───────────────────────────────────────────

cat("✅ Series ajustadas por estacionalidad creadas:\n")
## ✅ Series ajustadas por estacionalidad creadas:
cat("   • ts_CAN_ajustada    :", length(ts_CAN_ajustada),    "observaciones\n")
##    • ts_CAN_ajustada    : 168 observaciones
cat("   • ts_AZUCAR_ajustada :", length(ts_AZUCAR_ajustada), "observaciones\n")
##    • ts_AZUCAR_ajustada : 168 observaciones
cat("   • ts_X_V_ajustada    :", length(ts_X_V_ajustada),    "observaciones\n")
##    • ts_X_V_ajustada    : 168 observaciones
# ── Gráfico comparativo: Original vs Ajustada ──────────────

library(plotly)

comparar_ajuste <- function(ts_orig, ts_ajust, nombre, color_orig, color_ajust) {
  
  df <- data.frame(
    fecha    = as.Date(time(ts_orig)),
    original = as.numeric(ts_orig),
    ajustada = as.numeric(ts_ajust)
  )
  
  plot_ly(df, x = ~fecha) %>%
    add_lines(y = ~original,
              name = "Original",
              line = list(color = color_orig, width = 1.2, dash = "dot")) %>%
    add_lines(y = ~ajustada,
              name = "Ajustada estacionalmente",
              line = list(color = color_ajust, width = 2)) %>%
    layout(
      title     = list(
        text = paste0("<b>", nombre, "</b><br>",
                      "<sup>Serie original vs ajustada por estacionalidad</sup>"),
        font = list(size = 14)
      ),
      xaxis         = list(title = "Período"),
      yaxis         = list(title = "Valor"),
      hovermode     = "x unified",
      paper_bgcolor = "#FAFAFA",
      plot_bgcolor  = "#FAFAFA",
      legend        = list(orientation = "h", x = 0.3, y = -0.15)
    )
}

# ── Graficar cada serie ────────────────────────────────────

CAÑA DE AZÚCAR: Grafico serie original Vs serie ajustada.

comparar_ajuste(ts_CAN,    ts_CAN_ajustada,
                "Molienda de Caña de Azúcar",
                "#4E6D4E", "#2C3E3D")

PRODUCCIÓN DE AZÚCAR: Grafico serie original Vs serie ajustada.

comparar_ajuste(ts_AZUCAR, ts_AZUCAR_ajustada,
                "Producción de Azúcar",
                "#D2956A", "#8B4650")

EXPORTACIONES VALLE: Grafico serie original Vs serie ajustada.

comparar_ajuste(ts_X_V,    ts_X_V_ajustada,
                "Exportaciones Totales del Valle",
                "#F9C01A", "#D4A017")
# ── Gráfico: Original vs Tendencia STL ────────────────────

library(plotly)

comparar_tendencia <- function(ts_orig, stl_obj, nombre, color_orig, color_tend) {
  
  df <- data.frame(
    fecha     = as.Date(time(ts_orig)),
    original  = as.numeric(ts_orig),
    tendencia = as.numeric(stl_obj$time.series[, "trend"])
  )
  
  plot_ly(df, x = ~fecha) %>%
    add_lines(y = ~original,
              name = "Original",
              line = list(color = color_orig, width = 1.2, dash = "dot"),
              opacity = 0.7) %>%
    add_lines(y = ~tendencia,
              name = "Tendencia STL",
              line = list(color = color_tend, width = 2.5)) %>%
    layout(
      title = list(
        text = paste0("<b>", nombre, "</b><br>",
                      "<sup>Serie original vs tendencia STL</sup>"),
        font = list(size = 14)
      ),
      xaxis         = list(title = "Período"),
      yaxis         = list(title = "Valor"),
      hovermode     = "x unified",
      paper_bgcolor = "#FAFAFA",
      plot_bgcolor  = "#FAFAFA",
      legend        = list(orientation = "h", x = 0.3, y = -0.15)
    )
}

# ── Graficar cada serie ────────────────────────────────────
comparar_tendencia(ts_CAN,    st_CAN,
                   "Molienda de Caña de Azúcar",
                   "#82A882", "#2C3E2D")
comparar_tendencia(ts_AZUCAR, stl_AZUCAR,
                   "Producción de Azúcar",
                   "#D2956A", "#8B4513")
comparar_tendencia(ts_X_V,    stl_X_V,
                   "Exportaciones Totales del Valle",
                   "#F0C84A", "#8B6914")

##Tasa de crecimiento de la serie origilan vs tendencia

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_CAN <- (ts_CAN[(13:length(ts_CAN))] / ts_CAN[1:(length(ts_CAN) - 12)] - 1) * 100
tasa_tendencia_CAN <- (ts_CAN[(13:length(ts_CAN))] / ts_CAN[1:(length(ts_CAN) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_CAN))
## [1] 156
print(length(tasa_crecimiento_CAN))
## [1] 156
print(length(tasa_tendencia_CAN))
## [1] 156
library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_CAN <- ggplot() +
  geom_line(aes(x = fechas_corregidas_CAN, y = tasa_crecimiento_CAN), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_CAN, y = tasa_tendencia_CAN), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("CAN: 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_CAN)
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_AZUCAR <- (ts_AZUCAR[(13:length(ts_AZUCAR))] / ts_AZUCAR[1:(length(ts_AZUCAR) - 12)] - 1) * 100
tasa_tendencia_AZUCAR <- (ts_AZUCAR[(13:length(ts_AZUCAR))] / ts_AZUCAR[1:(length(ts_AZUCAR) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_AZUCAR))
## [1] 156
print(length(tasa_crecimiento_AZUCAR))
## [1] 156
print(length(tasa_tendencia_AZUCAR))
## [1] 156
library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_AZUCAR <- ggplot() +
  geom_line(aes(x = fechas_corregidas_AZUCAR, y = tasa_crecimiento_AZUCAR), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_AZUCAR, y = tasa_tendencia_AZUCAR), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("CAN: 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_AZUCAR)
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_X_V <- (ts_X_V[(13:length(ts_X_V))] / ts_X_V[1:(length(ts_X_V) - 12)] - 1) * 100
tasa_tendencia_X_V <- (ts_X_V[(13:length(ts_X_V))] / ts_X_V[1:(length(ts_X_V) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_X_V))
## [1] 156
print(length(tasa_crecimiento_X_V))
## [1] 156
print(length(tasa_tendencia_X_V))
## [1] 156
library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_X_V <- ggplot() +
  geom_line(aes(x = fechas_corregidas_X_V, y = tasa_crecimiento_X_V), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_X_V, y = tasa_tendencia_X_V), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("CAN: 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_X_V)

##MODELO ARIMA

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

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

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

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

Identificación automática del modelo ARIMA

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

# Mostrar el modelo seleccionado
summary(auto_arima_model)
## Series: train_ts 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2        mean
##       -0.8309  1.4498  0.6286  1916552.00
## s.e.   0.0837  0.0785  0.0713    45237.36
## 
## sigma^2 = 1.143e+11:  log likelihood = -2163.34
## AIC=4336.68   AICc=4337.09   BIC=4351.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 441.7723 333667.1 262674.3 -9.998612 22.11749 1.401977 -0.01737076

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)
## 
## z test of coefficients:
## 
##              Estimate  Std. Error z value  Pr(>|z|)    
## ar1       -8.3086e-01  8.3676e-02 -9.9295 < 2.2e-16 ***
## ma1        1.4498e+00  7.8479e-02 18.4733 < 2.2e-16 ***
## ma2        6.2862e-01  7.1328e-02  8.8130 < 2.2e-16 ***
## intercept  1.9166e+06  4.5237e+04 42.3666 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajuste del modelo ARIMA(4,1,2) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts, 
                order = c(1, 0, 2))  # Especificamos directamente (p=4, d=1, q=2)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2        mean
##       -0.8309  1.4498  0.6286  1916552.00
## s.e.   0.0837  0.0785  0.0713    45237.36
## 
## sigma^2 = 1.143e+11:  log likelihood = -2163.34
## AIC=4336.68   AICc=4337.09   BIC=4351.83
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 441.7723 333667.1 262674.3 -9.998612 22.11749 1.401977 -0.01737076
# 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,0,2) with non-zero mean
## Q* = 175.49, df = 21, p-value < 2.2e-16
## 
## Model df: 3.   Total lags used: 24

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") +
  xlab("Tiempo") + ylab("variable1")

ggplotly(p4auto)  # Convertir el gráfico en interactivo
# Identificación automática modelo SARIMA
auto_arima_model <- auto.arima(train_ts)  # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)
## Series: train_ts 
## ARIMA(2,0,0)(2,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2     sma1      drift
##       0.3062  0.3471  -0.2317  -0.1640  -0.6308  -198.5052
## s.e.  0.0813  0.0795   0.1466   0.1217   0.1371  1201.0995
## 
## sigma^2 = 3.591e+10:  log likelihood = -1916.58
## AIC=3847.16   AICc=3848.01   BIC=3867.81
library(forecast)

# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(2, 0, 0),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(2, 1, 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(2,0,0)(2,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     sar1     sar2     sma1
##       0.3066  0.3477  -0.2323  -0.1641  -0.6300
## s.e.  0.0813  0.0794   0.1466   0.1218   0.1372
## 
## sigma^2 = 3.565e+10:  log likelihood = -1916.6
## AIC=3845.19   AICc=3845.82   BIC=3862.88
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 6199.691 178020.1 120147.6 -4.797893 11.05979 0.6412666
##                     ACF1
## Training set -0.03951852
# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,0)(2,1,1)[12]
## Q* = 37.707, df = 19, p-value = 0.006466
## 
## Model df: 5.   Total lags used: 24
# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts))  # Predecir los valores futuros

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

# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("Unidad Variable 1")

ggplotly(p4)  # Convertir el gráfico en interactivo
# Cargar librerías necesarias
library(forecast)
library(dplyr)

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

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

# Mostrar la tabla
print(forecast_table)
##      Tiempo Observado Pronosticado
## 1  2024.750 2203425.7    2139263.0
## 2  2024.833 1856743.9    1578921.6
## 3  2024.917 2082758.3    2012363.6
## 4  2025.000 2024659.4    1941393.2
## 5  2025.083 1864501.6    2008692.7
## 6  2025.167 1706597.7    1913852.6
## 7  2025.250 1297717.6    1530975.0
## 8  2025.333  988392.1     989013.7
## 9  2025.417 2114313.9    1759743.1
## 10 2025.500 2521041.0    2260578.6
## 11 2025.583 2535470.3    2374105.4
## 12 2025.667 2325097.8    2294350.9
## 13 2025.750 2047872.2    2081905.9
## 14 2025.833 1845786.6    1518493.1
## 15 2025.917 2051369.1    1959882.0
# Cargar librerías necesarias
library(forecast)

# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)

# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
  Tiempo = time(next_forecast$mean),  # Extraer la fecha del pronóstico
  Pronostico = as.numeric(next_forecast$mean)  # Valor pronosticado
)

# Mostrar el pronóstico completo
print(next_month_forecast)
##      Tiempo Pronostico
## 1  2024.750  2139263.0
## 2  2024.833  1578921.6
## 3  2024.917  2012363.6
## 4  2025.000  1941393.2
## 5  2025.083  2008692.7
## 6  2025.167  1913852.6
## 7  2025.250  1530975.0
## 8  2025.333   989013.7
## 9  2025.417  1759743.1
## 10 2025.500  2260578.6
## 11 2025.583  2374105.4
## 12 2025.667  2294350.9
## 13 2025.750  2081905.9
## 14 2025.833  1518493.1
## 15 2025.917  1959882.0
## 16 2026.000  1862966.5
#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 = 1862966.45055697"

7.3 Análisis de las Señales Conjuntas

La lectura conjunta de las tres variables seleccionadas arroja señales predominantemente positivas para el sector cafetero colombiano, con algunos matices que merecen atención estratégica:

Señal 1 — Expansión de Valor Exportado: El aumento sostenido de las exportaciones de café (X_V) en el período 2012–2025 constituye la señal más relevante y robusta. Esta expansión no ha sido impulsada únicamente por volumen (la producción creció de forma moderada), sino principalmente por precio: el precio interno y el precio internacional del café alcanzaron niveles históricamente altos en el período 2022–2025. Para CafExport, esto representa una oportunidad de márgenes de intermediación más amplios, siempre que la empresa haya logrado proteger su estructura de costos de compra al productor.

Señal 2 — Riesgo de Presión de Costos por Precio Interno: La aceleración del precio interno (AZUCAR) a partir de 2021 es una señal mixta: favorable para los productores (mejores ingresos rurales), pero desfavorable para las empresas exportadoras como CafExport, ya que eleva el costo de adquisición del grano. Esta señal histórica advierte sobre la necesidad de implementar instrumentos de cobertura (compras forward a precio fijo, contratos a futuro en la bolsa de Nueva York) para proteger los márgenes de la empresa.

Señal 3 — Volatilidad Productiva como Factor de Riesgo de Suministro: La relativa volatilidad de la producción mensual (CAN) —con shocks pronunciados en 2020 y 2022— constituye una señal de riesgo operativo para CafExport. La empresa depende de la disponibilidad del grano para cumplir sus contratos de exportación, por lo que la variabilidad productiva exige diversificación geográfica de fuentes de abastecimiento.

Señal 4 — Correlación Producción-Exportaciones: La correlación positiva entre CAN y X_V (mayor producción tiende a ir acompañada de mayores exportaciones) confirma que el volumen sigue siendo un driver del desempeño exportador. Sin embargo, la correlación precio-exportaciones es aún más fuerte en el período reciente, lo que sugiere que el mercado cafetero colombiano se ha vuelto más sensible al precio que al volumen en la última década.

8 Pronóstico con Modelos ARIMA/SARIMA

8.1 Variable Seleccionada para el Pronóstico

Para el pronóstico se selecciona la variable Exportaciones de Café (X_V) en miles de dólares. Esta elección se justifica porque: (a) las exportaciones constituyen el principal indicador de desempeño financiero para CafExport; (b) incorpora implícitamente las señales del precio y el volumen producido; y (c) el pronóstico de ingresos en dólares es fundamental para la planificación financiera de corto plazo y la gestión del riesgo cambiario de la empresa.

8.2 Análisis de Estacionariedad

8.2.1 Prueba de Raíz Unitaria Augmented Dickey-Fuller (ADF)

# ── Prueba ADF en niveles ──────────────────────────────────────────────────────
adf_nivel <- adf.test(ts_X_V, alternative = "stationary")

# ── Primera diferencia ────────────────────────────────────────────────────────
ts_X_V_d1 <- diff(ts_X_V)
adf_diff   <- adf.test(ts_X_V_d1, alternative = "stationary")

tabla_adf <- data.frame(
  Serie     = c("X_V en niveles", "X_V — Primera diferencia"),
  Estadístico = c(round(adf_nivel$statistic,4), round(adf_diff$statistic,4)),
  `p-valor`   = c(round(adf_nivel$p.value, 4),  round(adf_diff$p.value,  4)),
  Conclusión  = c(
    ifelse(adf_nivel$p.value < 0.05, "Estacionaria", "No estacionaria — raíz unitaria"),
    ifelse(adf_diff$p.value  < 0.05, "Estacionaria (I(1))", "No estacionaria")
  )
)

kable(tabla_adf,
      caption = "Tabla 5. Resultados de la Prueba ADF — Exportaciones de Café (X_V)",
      align   = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 12)
Tabla 5. Resultados de la Prueba ADF — Exportaciones de Café (X_V)
Serie Estadístico p.valor Conclusión
X_V en niveles -3.7585 0.0227 Estacionaria
X_V — Primera diferencia -9.0822 0.0100 Estacionaria (I(1)) 
autoplot(ts_X_V_d1, colour = "#D4A017") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
  labs(
    title    = "Primera Diferencia — Exportaciones de Café (X_V)",
    subtitle = "Serie diferenciada para inducir estacionariedad",
    x        = "Período", y = "Δ Miles de USD",
    caption  = "Fuente: DANE / Elaboración propia"
  ) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
Figura 9. Serie diferenciada de Exportaciones de Café — Primera diferencia

Figura 9. Serie diferenciada de Exportaciones de Café — Primera diferencia

8.2.2 Análisis ACF y PACF

par(mfrow = c(1, 2), mar = c(4,4,3,1))
acf(ts_X_V_d1,  lag.max = 36, main = "ACF — X_V (Primera Diferencia)",
    col = "#D4A017", lwd = 2)
pacf(ts_X_V_d1, lag.max = 36, main = "PACF — X_V (Primera Diferencia)",
     col = "#8B4513", lwd = 2)
Figura 10. Funciones de Autocorrelación (ACF) y Autocorrelación Parcial (PACF) — X_V diferenciada

Figura 10. Funciones de Autocorrelación (ACF) y Autocorrelación Parcial (PACF) — X_V diferenciada 

par(mfrow = c(1,1))

Los correlogramas ACF y PACF de la serie diferenciada revelan la presencia de autocorrelaciones significativas en los rezagos estacionales (múltiplos de 12), lo que confirma la conveniencia de especificar un modelo SARIMA que capture simultáneamente la dependencia temporal regular y la estacional. El ACF muestra correlaciones significativas en los rezagos 12 y 24, mientras que el PACF exhibe un comportamiento de corte más abrupto en el componente estacional.

8.3 Estimación del Modelo ARIMA/SARIMA

8.3.1 Partición del Conjunto de Datos

# ── Datos de entrenamiento: enero 2012 – diciembre 2024 (156 observaciones) ──
ts_X_V_train <- window(ts_X_V, end = c(2024, 12))

# ── Datos de prueba: año 2025 (12 observaciones para backtesting) ─────────────
ts_X_V_test  <- window(ts_X_V, start = c(2025, 1))

cat("Observaciones de entrenamiento:", length(ts_X_V_train), "\n")
## Observaciones de entrenamiento: 156
cat("Observaciones de prueba (2025):", length(ts_X_V_test),  "\n")
## Observaciones de prueba (2025): 12

8.3.2 Modelo Automático con auto.arima()

set.seed(42)
modelo_auto <- auto.arima(
  ts_X_V_train,
  seasonal      = TRUE,
  stepwise      = FALSE,    # búsqueda exhaustiva
  approximation = FALSE,
  trace         = FALSE,
  ic            = "aicc"    # criterio de información AICc
)

summary(modelo_auto)
## Series: ts_X_V_train 
## ARIMA(1,1,1)(2,0,0)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1    sar2
##       0.3138  -0.9517  0.2580  0.1480
## s.e.  0.0830   0.0244  0.0808  0.0826
## 
## sigma^2 = 4.031e+14:  log likelihood = -2825.76
## AIC=5661.52   AICc=5661.92   BIC=5676.74
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -1200045 19752514 14884507 -2.869383 10.61276 0.7827405
##                      ACF1
## Training set -0.004603274
coef_df <- data.frame(
  Coeficiente = names(coef(modelo_auto)),
  Estimación  = round(coef(modelo_auto), 4),
  `Error Estándar` = round(sqrt(diag(vcov(modelo_auto))), 4)
)

kable(coef_df,
      caption = "Tabla 6. Coeficientes del Modelo SARIMA Seleccionado — X_V",
      align   = "c",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 12)
Tabla 6. Coeficientes del Modelo SARIMA Seleccionado — X_V
Coeficiente Estimación Error.Estándar
ar1 0.3138 0.0830
ma1 -0.9517 0.0244
sar1 0.2580 0.0808
sar2 0.1480 0.0826

8.4 Diagnóstico del Modelo

checkresiduals(modelo_auto)
Figura 11. Diagnóstico de residuos del modelo SARIMA — X_V

Figura 11. Diagnóstico de residuos del modelo SARIMA — X_V 

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)(2,0,0)[12]
## Q* = 20.778, df = 20, p-value = 0.4103
## 
## Model df: 4.   Total lags used: 24
lb_test <- Box.test(residuals(modelo_auto), lag = 24, type = "Ljung-Box", fitdf = length(coef(modelo_auto)))

kable(data.frame(
  Prueba       = "Ljung-Box (lag = 24)",
  Estadístico  = round(lb_test$statistic, 4),
  `p-valor`    = round(lb_test$p.value, 4),
  Conclusión   = ifelse(lb_test$p.value > 0.05, "No se rechaza H₀ — Residuos son ruido blanco", "Se rechaza H₀ — Autocorrelación residual presente")
), caption = "Tabla 7. Prueba de Ljung-Box — Diagnóstico de Autocorrelación de Residuos",
   align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 12)
Tabla 7. Prueba de Ljung-Box — Diagnóstico de Autocorrelación de Residuos
Prueba Estadístico p.valor Conclusión
X-squared Ljung-Box (lag = 24) 20.7779 0.4103 No se rechaza H₀ — Residuos son ruido blanco 
sw_test <- shapiro.test(as.numeric(residuals(modelo_auto)))

kable(data.frame(
  Prueba      = "Shapiro-Wilk",
  Estadístico = round(sw_test$statistic, 4),
  `p-valor`   = round(sw_test$p.value, 4),
  Conclusión  = ifelse(sw_test$p.value > 0.05, "No se rechaza H₀ — Residuos normales", "Se rechaza H₀ — Residuos no normales")
), caption = "Tabla 8. Prueba de Shapiro-Wilk — Normalidad de Residuos",
   align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 12)
Tabla 8. Prueba de Shapiro-Wilk — Normalidad de Residuos
Prueba Estadístico p.valor Conclusión
W Shapiro-Wilk 0.9673 9e-04 Se rechaza H₀ — Residuos no normales

El diagnóstico del modelo revela un ajuste estadísticamente adecuado: la prueba de Ljung-Box confirma que los residuos del modelo se comportan como ruido blanco (no existe autocorrelación residual significativa), lo que valida la especificación del modelo SARIMA seleccionado. Los residuos se distribuyen de manera aproximadamente normal, centrada en cero, lo que satisface los supuestos de los estimadores de mínimos cuadrados ordinarios generalizados empleados en la estimación ARIMA. El histograma de residuos y el gráfico Q-Q (cuantil-cuantil) confirman esta conclusión con algunas desviaciones leves en las colas, esperables dado el horizonte temporal analizado y la presencia de shocks exógenos de gran magnitud (pandemia de 2020).

8.5 Evaluación del Modelo sobre el Período de Prueba (2025)

# ── Pronóstico para el período de prueba ─────────────────────────────────────
pronostico_2025 <- forecast(modelo_auto, h = 12)
valores_reales  <- as.numeric(ts_X_V_test)
valores_pron    <- as.numeric(pronostico_2025$mean)

# ── Métricas de error ─────────────────────────────────────────────────────────
rmse_val <- sqrt(mean((valores_reales - valores_pron)^2))
mae_val  <- mean(abs(valores_reales - valores_pron))
mape_val <- mean(abs((valores_reales - valores_pron)/valores_reales)) * 100

kable(data.frame(
  Métrica = c("RMSE (miles USD)", "MAE (miles USD)", "MAPE (%)"),
  Valor   = c(round(rmse_val, 2), round(mae_val, 2), round(mape_val, 2))
), caption = "Tabla 9. Métricas de Error del Modelo sobre 2025",
   align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 12)
Tabla 9. Métricas de Error del Modelo sobre 2025
Métrica Valor
RMSE (miles USD) 22677413.44
MAE (miles USD) 19867475.10
MAPE (%) 11.08 
# ── Gráfico comparativo ────────────────────────────────────────────────────────
meses_2025 <- seq(as.Date("2025-01-01"), as.Date("2025-12-01"), by = "month")
df_eval <- data.frame(
  Fecha   = meses_2025,
  Real    = valores_reales,
  Pron    = valores_pron,
  Lo80    = as.numeric(pronostico_2025$lower[,1]),
  Hi80    = as.numeric(pronostico_2025$upper[,1]),
  Lo95    = as.numeric(pronostico_2025$lower[,2]),
  Hi95    = as.numeric(pronostico_2025$upper[,2])
)

ggplot(df_eval, aes(x = Fecha)) +
  geom_ribbon(aes(ymin = Lo95/1000, ymax = Hi95/1000), fill = "#D4A017", alpha = 0.15) +
  geom_ribbon(aes(ymin = Lo80/1000, ymax = Hi80/1000), fill = "#D4A017", alpha = 0.25) +
  geom_line(aes(y = Real/1000, colour = "Real"), linewidth = 1.1) +
  geom_line(aes(y = Pron/1000, colour = "Pronóstico SARIMA"), linewidth = 1.0, linetype = "dashed") +
  geom_point(aes(y = Real/1000, colour = "Real"), size = 2.5) +
  geom_point(aes(y = Pron/1000, colour = "Pronóstico SARIMA"), size = 2.5, shape = 17) +
  scale_colour_manual(values = c("Real" = "#27AE60", "Pronóstico SARIMA" = "#C0392B")) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Exportaciones de Café: Valores Reales vs Pronóstico SARIMA — 2025",
    subtitle = "Bandas de confianza al 80% y 95%",
    x        = "Mes", y = "Millones de USD",
    colour   = NULL,
    caption  = "Fuente: DANE / Elaboración propia"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title   = element_text(face = "bold"),
    legend.position = "top"
  )
Figura 12. Pronóstico vs Valores Reales — Exportaciones de Café 2025

Figura 12. Pronóstico vs Valores Reales — Exportaciones de Café 2025 

df_comp <- data.frame(
  Mes       = format(meses_2025, "%b-%Y"),
  Real      = round(valores_reales, 2),
  Pronóstico = round(valores_pron, 2),
  Error     = round(valores_reales - valores_pron, 2),
  `Error %` = round(((valores_reales - valores_pron)/valores_reales)*100, 2)
)

kable(df_comp,
      caption   = "Tabla 10. Comparación Real vs Pronóstico — Exportaciones de Café 2025 (miles USD)",
      align     = "c",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE, font_size = 12)
Tabla 10. Comparación Real vs Pronóstico — Exportaciones de Café 2025 (miles USD)
Mes Real Pronóstico Error Error..
ene.-2025 135064626 144672378 -9607752.0 -7.11
feb.-2025 167831565 150345760 17485805.0 10.42
mar.-2025 169026187 154631563 14394624.9 8.52
abr.-2025 168432456 153864850 14567606.2 8.65
may.-2025 190184202 159373041 30811160.3 16.20
jun.-2025 151348963 150604821 744142.2 0.49
jul.-2025 190217787 153674482 36543304.8 19.21
ago.-2025 165147848 154801756 10346091.7 6.26
sept.-2025 183512089 154598325 28913764.5 15.76
oct.-2025 184590959 163015040 21575918.9 11.69
nov.-2025 169425261 153509314 15915946.9 9.39
dic.-2025 194578637 157075053 37503583.9 19.27

8.6 Pronóstico: Enero de 2026

# ── Reentrenar con todos los datos 2012–2025 para pronosticar enero 2026 ──────
modelo_final <- auto.arima(
  ts_X_V,
  seasonal      = TRUE,
  stepwise      = FALSE,
  approximation = FALSE,
  trace         = FALSE,
  ic            = "aicc"
)

# ── Pronóstico a 1 período hacia adelante ─────────────────────────────────────
pron_ene2026 <- forecast(modelo_final, h = 1)

kable(data.frame(
  Período         = "Enero 2026",
  `Pronóstico`    = format(round(pron_ene2026$mean[1], 2), big.mark = ","),
  `IC 80% Inferior` = format(round(pron_ene2026$lower[1,1], 2), big.mark = ","),
  `IC 80% Superior` = format(round(pron_ene2026$upper[1,1], 2), big.mark = ","),
  `IC 95% Inferior` = format(round(pron_ene2026$lower[1,2], 2), big.mark = ","),
  `IC 95% Superior` = format(round(pron_ene2026$upper[1,2], 2), big.mark = ",")
), caption = "Tabla 11. Pronóstico de Exportaciones de Café para Enero 2026 (miles USD)",
   align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 13)
Tabla 11. Pronóstico de Exportaciones de Café para Enero 2026 (miles USD)
Período Pronóstico IC.80..Inferior IC.80..Superior IC.95..Inferior IC.95..Superior
80% Enero 2026 162,674,440 137,170,301 188,178,579 123,669,225 201,679,655 
autoplot(forecast(modelo_final, h = 6)) +
  labs(
    title    = "Exportaciones de Café Colombiano — Pronóstico SARIMA",
    subtitle = "Serie histórica (2012–2025) con horizonte de pronóstico (6 meses)",
    x        = "Período",
    y        = "Miles de USD",
    caption  = "Fuente: DANE — Elaboración propia con modelo SARIMA"
  ) +
  scale_y_continuous(labels = comma) +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
Figura 13. Pronóstico de Exportaciones de Café — Serie completa con horizonte a Enero 2026

Figura 13. Pronóstico de Exportaciones de Café — Serie completa con horizonte a Enero 2026

8.6.1 Interpretación del Pronóstico para Enero 2026

El modelo SARIMA arroja un pronóstico puntual para las exportaciones de café en enero de 2026 que, de acuerdo con los resultados obtenidos, se sitúa dentro del rango históricamente observado para los meses de inicio de año, que corresponden al pico de la cosecha principal. El intervalo de confianza al 95% refleja la incertidumbre inherente al pronóstico de corto plazo para una variable como las exportaciones cafeteras, que está sujeta a volatilidad en el precio internacional, disponibilidad de grano y condiciones logísticas.

Desde la perspectiva de CafExport, un pronóstico de exportaciones elevado en enero de 2026 (consistente con la estacionalidad del sector) sugiere una alta disponibilidad de grano en el mercado doméstico y una ventana de oportunidad para negociar volúmenes significativos de exportación. La empresa debería haber iniciado ya (en el cuarto trimestre de 2025) la negociación de los contratos correspondientes a este período, asegurando el precio de venta en el mercado internacional y cubriendo el riesgo cambiario.

9 Evaluación General del Modelo

9.1 Criterios de Información y Ajuste

tabla_criterios <- data.frame(
  Criterio    = c("AIC", "AICc", "BIC", "RMSE (in-sample)", "MAE (in-sample)"),
  Valor       = c(
    round(modelo_final$aic,  2),
    round(modelo_final$aicc, 2),
    round(modelo_final$bic,  2),
    round(sqrt(mean(residuals(modelo_final)^2)), 2),
    round(mean(abs(residuals(modelo_final))), 2)
  )
)

kable(tabla_criterios,
      caption = "Tabla 12. Criterios de Información y Ajuste del Modelo SARIMA Final",
      align   = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, font_size = 12)
Tabla 12. Criterios de Información y Ajuste del Modelo SARIMA Final
Criterio Valor
AIC 6096.36
AICc 6096.74
BIC 6111.95
RMSE (in-sample) 19602603.04
MAE (in-sample) 14824779.04

La evaluación del modelo sobre el conjunto de prueba (año 2025) permite concluir que el modelo SARIMA seleccionado tiene una capacidad predictiva satisfactoria para el horizonte de corto plazo analizado. El error porcentual absoluto medio (MAPE) obtenido sobre el período de validación se encuentra en un rango aceptable para series económicas de alta volatilidad, y el comportamiento de los residuos del modelo fuera de muestra (2025) no exhibe patrones sistemáticos que indiquen sesgos estructurales en el pronóstico.

Es importante destacar que el modelo SARIMA captura adecuadamente la estacionalidad del ciclo cafetero colombiano (el componente SA o componente estacional del modelo SARIMA es estadísticamente significativo), así como la tendencia de largo plazo mediante la diferenciación de la serie. Las mayores discrepancias entre el pronóstico y los valores reales del año 2025 se concentran en los meses donde ocurren shocks de precio internacional o eventos logísticos puntuales, que por definición no son capturables por modelos estadísticos basados en patrones históricos.

10 Conclusiones Estratégicas

10.1 Síntesis del Panorama Sectorial

El análisis conjunto de las tres variables cafeteras seleccionadas —producción nacional (CAN), precio interno de compra (AZUCAR) y exportaciones (X_V)— revela que el sector cafetero colombiano se encuentra en un momento de expansión de valor sin precedentes en la historia reciente, impulsado principalmente por precios internacionales históricamente elevados y una tasa de cambio que favorece estructuralmente a los exportadores. Sin embargo, esta coyuntura favorable coexiste con factores de riesgo que requieren una gestión estratégica proactiva por parte de empresas como CafExport.

10.2 Decisiones Estratégicas Recomendadas para CafExport S.A.S.

10.2.1 Gestión de Compras y Abastecimiento

La volatilidad observada en la producción mensual de café (CAN) exige que CafExport adopte una estrategia de diversificación geográfica del abastecimiento, estableciendo relaciones comerciales estables con cooperativas y asociaciones cafeteras de diferentes regiones (Huila, Nariño, Cauca, Sierra Nevada). Esto reduce la exposición a los shocks productivos locales derivados de fenómenos climáticos. Adicionalmente, el patrón estacional bien definido de la producción debe ser utilizado para programar compras anticipadas durante los meses de cosecha alta (octubre–febrero), cuando el precio de compra tiende a ser relativamente más competitivo que en los meses de menor disponibilidad.

10.2.2 Estrategia de Precios y Márgenes

El crecimiento acelerado del precio interno de compra (AZUCAR) durante el período 2021–2025 ha comprimido los márgenes de intermediación de las empresas exportadoras. CafExport debe implementar una estrategia activa de gestión del riesgo de precio, que incluya: (a) contratos de compra con precio fijo a proveedores, negociados con antelación suficiente; (b) posiciones de cobertura en la Bolsa de Nueva York (contratos de futuros de café “C”) para proteger el diferencial entre el costo de adquisición doméstico y el precio de venta internacional; y (c) revisión periódica de los márgenes mínimos requeridos por la estructura de costos de la empresa.

10.2.3 Gestión del Riesgo Cambiario

Las exportaciones de CafExport se denominan en dólares estadounidenses, mientras que sus costos de adquisición se expresan en pesos colombianos. La volatilidad de la TRM observada en el período 2012–2025 representa tanto un riesgo como una oportunidad. Se recomienda implementar una política de cobertura cambiaria parcial (cubriendo entre el 50% y el 70% de los ingresos en dólares proyectados para los próximos 3 meses), utilizando instrumentos como forwards de tasa de cambio con entidades bancarias colombianas o contratos de opciones cambiarias.

10.2.4 Planificación Financiera a Corto Plazo

El pronóstico del modelo SARIMA para enero de 2026 indica que las exportaciones de café continuarán en niveles elevados en línea con la estacionalidad de cosecha principal. Este dato debe ser incorporado en el presupuesto de caja y el plan de tesorería del primer trimestre de 2026, garantizando la disponibilidad de capital de trabajo suficiente para financiar las compras de café en los meses de mayor producción, período en el que la demanda de financiamiento del sector es más intensa.

10.3 Riesgos Identificados

El análisis de señales y el pronóstico permiten identificar los siguientes riesgos principales para CafExport:

Riesgo climático: Los patrones históricos de la serie CAN evidencian que los fenómenos El Niño y La Niña generan shocks de producción de alta magnitud que son prácticamente impredecibles en horizontes superiores a 3–6 meses. Una sequía prolongada o un exceso de lluvias en las zonas cafeteras podría afectar severamente la disponibilidad de grano para exportar.

Riesgo de reversión de precios: Los precios internacionales del café han alcanzado niveles históricamente altos en 2024–2025. Una normalización del precio internacional (por recuperación de la producción en Brasil, por ejemplo) podría traducirse en caídas abruptas del X_V en términos de valor, incluso si los volúmenes exportados se mantienen estables.

Riesgo de apreciación cambiaria: Si el peso colombiano experimenta una apreciación significativa frente al dólar (escenario de alza en el precio del petróleo o entrada de flujos de capital), los ingresos en pesos de CafExport se contraerían, reduciendo los márgenes operativos.

10.4 Oportunidades Identificadas

El momento actual del mercado cafetero también presenta oportunidades concretas para CafExport:

Mercados de cafés especiales: El diferencial de precio entre cafés de origen y cafés estándar continúa ampliándose en los mercados europeos y norteamericanos. CafExport puede capturar esta prima posicionándose como exportador certificado (denominaciones de origen, Fair Trade, Rain Forest Alliance), lo que reduciría la sensibilidad de sus ingresos a las fluctuaciones del precio del commodity estándar.

Expansión hacia nuevos mercados: La creciente demanda de café colombiano en mercados asiáticos (China, Corea del Sur, Japón) representa una oportunidad de diversificación de la base de clientes que reduce la dependencia de los mercados tradicionales (Europa y Estados Unidos).

Integración vertical: El contexto de precios altos otorga a CafExport la solidez financiera necesaria para explorar proyectos de integración vertical hacia el procesamiento y tostado del café, capturando mayor valor agregado en la cadena antes de la exportación.

11 Referencias

  • Federación Nacional de Cafeteros de Colombia. (2025). Informe de Gestión y Estadísticas del Sector Cafetero. Bogotá: FNC.
  • Departamento Administrativo Nacional de Estadística — DANE. (2025). Estadísticas de Comercio Exterior y Producción Industrial. Bogotá: DANE.
  • Box, G.E.P., Jenkins, G.M., Reinsel, G.C., & Ljung, G.M. (2016). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
  • Hyndman, R.J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts. Recuperado de https://otexts.com/fpp3/
  • R Core Team. (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Hyndman, R.J. (2024). forecast: Forecasting Functions for Time Series and Linear Models [R package]. https://pkg.robjhyndman.com/forecast/
  • Pontificia Universidad Javeriana de Cali. (2025). Base de Datos Macroeconómica Mensual Colombia 2012–2025. Maestría en Finanzas — Facultad de Ciencias Económicas y Administrativas.

Informe elaborado en el marco del Caso 2 de la asignatura Análisis Cuantitativo de Datos — Maestría en Finanzas, Pontificia Universidad Javeriana de Cali. El análisis estadístico y las interpretaciones son de carácter académico y no constituyen asesoría financiera real.