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("Base Caso2.xlsx")
View(data_col)

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

Variable 1

# Convertir/declarar variable 1=IPC exportacion de cafe   en serie de tiempo mensual  
variable1_ts <- ts(data_col$IPC, start = c(2012, 1), frequency = 12)

Variable 2

# Convertir/declarar el PNCAFE produccion de cafe en serie de tiempo  mensual
variable2_ts <- ts(data_col$PNCAFE, start = c(2012, 1), frequency = 12)

Variable 3

# Convertir/declarar PECAFE precio externo de cafe col en serie de tiempo mensual
variable3_ts <- ts(data_col$PECAFE, start = c(2012, 1), frequency = 12)

#1. Introducción

El sector elegido es el agroexportador, con énfasis en la industria del café en Colombia, dado su relevancia económica, social y comercial. La empresa escogida es Cafexport S.A., dedicada a la producción y exportación de café.

Las tres variables seleccionadas son:

IPC: índice de precios al consumidor (contexto macroeconómico e inflación interna).

PNCAFE: indicador nacional relacionado con el café (precio o producción a nivel interno).

PECAFE: precio de exportación del café (variable clave para ingresos y competitividad externa).

Justificación: El IPC impacta los costos internos y el poder adquisitivo; PNCAFE refleja la dinámica de la cadena de suministro nacional; PECAFE define los ingresos en mercados internacionales. La combinación de las tres permite evaluar el desempeño y perspectivas de la empresa.

Metodologia

Se aplicó la descomposición STL (Seasonal-Trend decomposition using Loess) a cada variable, extrayendo los componentes de tendencia, estacionalidad y residuo. Posteriormente, se aplicaron pruebas de estacionariedad (ADF) y se estimaron modelos ARIMA automáticos (con y sin componente estacional, SARIMA) sobre la variable PECAFE, seleccionada por su impacto directo en ingresos por exportaciones. Finalmente, se generó un pronóstico para enero de 2025 y se evaluó su precisión mediante comparación con valores recientes.

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 = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Extracción señales variable 1

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

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

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

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

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

#analisis Tendencia: crecimiento sostenido y estable a lo largo del período (2012–2024), lo cual refleja el aumento generalizado y acumulativo de precios en la economía.

Estacionalidad: patrón repetitivo anual, con fluctuaciones relativamente pequeñas alrededor de la tendencia.

Residuo: variaciones de baja magnitud, sin choques extremos, lo que confirma que el IPC es una serie bastante estable.

Interpretación: el IPC muestra el proceso inflacionario estructural del país, lo que afecta los costos internos de producción y el poder adquisitivo de los consumidores.

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)

#analisis

Tendencia: crecimiento inicial fuerte hasta 2015, luego una fase de estancamiento/caída entre 2016–2020 y finalmente una recuperación clara en los últimos años.

Estacionalidad: muy marcada, con picos y valles anuales pronunciados, reflejando la estacionalidad de la producción y comercialización del café.

Residuo: alta volatilidad, lo que sugiere choques externos (clima, plagas, variaciones de oferta y demanda).

Interpretación: PNCAFE refleja la sensibilidad de la producción/precio nacional a factores coyunturales y estacionales, mostrando mayor inestabilidad que el IPC.

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)

#analisis Tendencia: muestra fases de caída (2012–2015), estancamiento (2016–2019) y recuperación significativa desde 2020, llegando a máximos en 2024.

Estacionalidad: presente y regular, aunque con menor amplitud que PNCAFE.

Residuo: choques fuertes en años específicos, lo que puede asociarse con crisis internacionales de precios de café o cambios en la demanda externa.

Interpretación: PECAFE depende fuertemente del mercado internacional y refleja los ingresos por exportaciones. La recuperación reciente coincide con mejores precios internacionales.

#comparacion

IPC vs PNCAFE/PECAFE: mientras el IPC muestra una trayectoria estable y creciente, PNCAFE y PECAFE son más volátiles y responden a factores sectoriales.

PNCAFE vs PECAFE: ambas están correlacionadas; los picos y valles del PNCAFE suelen coincidir con los de PECAFE, lo que confirma la conexión entre el mercado interno y el externo.

Tasa de crecimiento anual (trend):

IPC → estable, entre 3–6 % anual aprox.

PNCAFE → tasas muy irregulares, con periodos negativos (caídas) y otros de fuerte crecimiento.

PECAFE → negativas hasta 2015, luego positivas y actualmente en un ciclo alcista.

#Insights

Estabilidad macro vs volatilidad sectorial: el IPC sube de manera estable, mientras que PNCAFE y PECAFE son más inestables, reflejando que el sector cafetero está sujeto a riesgos propios (clima, mercado internacional).

Dependencia externa: PECAFE marca los ingresos por exportaciones y determina la rentabilidad de las empresas cafeteras, mostrando un repunte favorable desde 2020.

Correlación interna-externa: la fuerte similitud entre PNCAFE y PECAFE confirma que los precios nacionales siguen la dinámica internacional.

Oportunidad estratégica: el crecimiento reciente del PECAFE puede aprovecharse para asegurar contratos a precios altos, mientras que la estacionalidad y volatilidad sugieren la necesidad de coberturas y planificación financiera.

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 = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 1:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)

Gráfico serie original VS ajustada Variable 2

# Crear vector de fechas correctamente alineado con la serie
fechas_var2 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var2 <- ggplot() +
  geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 2:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 2") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)

Gráfico serie original VS ajustada Variable 3

# Crear vector de fechas correctamente alineado con la serie
fechas_var3 <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))

# Gráfico mejorado con fechas en el eje X
grafico_ajustada_var3 <- ggplot() +
  geom_line(aes(x = fechas_var3, y = variable3_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 3:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 3") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)

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" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 1: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

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" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 2: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 2") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

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" = "grey", "Tendencia" = "black")) +
  ggtitle("Variable 3: Serie Original vs Tendencia") +
  xlab("Tiempo") +
  ylab("Unidad de medida Variable 3") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X

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

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] 144
print(length(tasa_crecimiento_var1))
## [1] 144
print(length(tasa_tendencia_var1))
## [1] 144

*Gráfico variable original y tendencia variable 1: tasa de crecimiento anual**

library(ggplot2)
library(plotly)

# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable1: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

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] 144
print(length(tasa_crecimiento_var2))
## [1] 144
print(length(tasa_tendencia_var2))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable2: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

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] 144
print(length(tasa_crecimiento_var3))
## [1] 144
print(length(tasa_tendencia_var3))
## [1] 144
# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var3 <- ggplot() +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
  ggtitle("Variable3: Tasa de crecimiento anual % de la serie Original y la tendencia") +
  xlab("Tiempo") +
  ylab("% de Crecimiento Anual") +
  theme_minimal()

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

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 2024. Conjunto de prueba (test): Datos desde octubre de 2024 hasta diciembre de 2024.

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

# 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(variable1_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable1_ts, end = c(2024, 9))  # Entrenamiento hasta septiembre 2024
test_ts <- window(variable1_ts, start = c(2024, 10))  # Prueba inicia desde oct2024

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(2,2,1) 
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.8153  -0.2195  -0.9082
## s.e.  0.0949   0.0861   0.0622
## 
## sigma^2 = 0.0801:  log likelihood = -22.47
## AIC=52.94   AICc=53.21   BIC=65.01
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE       MASE
## Training set 0.01061054 0.2783565 0.2042477 0.01040633 0.1984684 0.03786293
##                      ACF1
## Training set -0.008880447

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|)    
## ar1  0.815306   0.094904   8.5908  < 2e-16 ***
## ar2 -0.219452   0.086113  -2.5484  0.01082 *  
## ma1 -0.908227   0.062167 -14.6094  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajuste del modelo ARIMA(2,2,1) 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(2, 2, 1))  # Especificamos directamente (p=2, d=2, q=1)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)
## Series: train_ts 
## ARIMA(2,2,1) 
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.8153  -0.2195  -0.9082
## s.e.  0.0949   0.0861   0.0622
## 
## sigma^2 = 0.0801:  log likelihood = -22.47
## AIC=52.94   AICc=53.21   BIC=65.01
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE       MASE
## Training set 0.01061054 0.2783565 0.2042477 0.01040633 0.1984684 0.03786293
##                      ACF1
## Training set -0.008880447

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(2,2,1)
## Q* = 58.905, df = 21, p-value = 1.868e-05
## 
## 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

Interpretación modelo automatico (2,2,1):El modelo automático (2,2,1 ) parece pronosticar mejor dentro de prueba. Hay un sobre ajuste, pero se capturan muy bien los puntos de quiebre. Es un modelo tentativo adecuado para pronpostico fuera de muestra o a futuro.

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 2024.750    143.83     144.5871
## 2 2024.833    144.22     145.2544
## 3 2024.917    144.88     145.9558

Análisis del modelo ARIMA (2,2,1)

Diagnóstico de residuos

Los residuos del modelo (gráfico 1) se distribuyen alrededor de cero, sin patrón claro → indican que el modelo captura bastante bien la estructura de la serie.

El histograma de los residuos se aproxima a una distribución normal, aunque con colas algo pesadas.

La función de autocorrelación (ACF) no muestra lags significativamente fuera de las bandas de confianza → evidencia de que no quedaron correlaciones importantes sin modelar. Conclusión: el modelo es estadísticamente adecuado.

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

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

# Cargar librerías necesarias
library(forecast)

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

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

# Mostrar el pronóstico completo
print(next_month_forecast_auto)
##     Tiempo Pronostico
## 1 2024.750   144.5871
## 2 2024.833   145.2544
## 3 2024.917   145.9558
## 4 2025.000   146.6630
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2025 = 146.662959097937"

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 dautomá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(1,1,2)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1      ma2     sma1
##       0.9537  -0.3101  -0.2039  -0.8243
## s.e.  0.0352   0.0895   0.0866   0.0866
## 
## sigma^2 = 0.05467:  log likelihood = 0.72
## AIC=8.57   AICc=9.01   BIC=23.27

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

# Cargar el paquete necesario
library(forecast)

# Ajustar el modelo SARIMA(1,1,2)(0,1,1)[12]  #Modelo identificado en el paso anterior
darima <- Arima(train_ts, 
                order = c(1, 1, 2),  # (p,d,q) -> (1,1,2)
                seasonal = list(order = c(0, 1, 1),  # (P,D,Q) -> (0,1,1)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)
## Series: train_ts 
## ARIMA(1,1,2)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1      ma2     sma1
##       0.9537  -0.3101  -0.2039  -0.8243
## s.e.  0.0352   0.0895   0.0866   0.0866
## 
## sigma^2 = 0.05467:  log likelihood = 0.72
## AIC=8.57   AICc=9.01   BIC=23.27
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE       MASE
## Training set 0.01550025 0.2204357 0.1555902 0.01633994 0.1502826 0.02884292
##                    ACF1
## Training set 0.01873309

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(1,1,2)(0,1,1)[12]
## Q* = 21.057, df = 20, p-value = 0.3938
## 
## Model df: 4.   Total lags used: 24

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") +
  xlab("Tiempo") + ylab("Unidad Variable 1")

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

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 2024.750    143.83     144.2030
## 2 2024.833    144.22     144.4865
## 3 2024.917    144.88     145.0330

Analisis

El diagnóstico del modelo ARIMA aplicado al IPC, la producción de café (PNCAFE) y el precio de exportación de café (PECAFE) muestra resultados en general satisfactorios. Los residuos presentan un comportamiento cercano a un ruido blanco, con distribución aproximadamente normal y sin autocorrelaciones significativas relevantes, lo que respalda la validez del modelo. Si bien existen ligeras desviaciones en las colas y algunos rezagos cercanos a los límites, estos no comprometen de manera significativa la robustez del pronóstico. En la comparación entre valores observados y pronosticados, el modelo logra reproducir adecuadamente la dirección de la tendencia. Tanto en la gráfica como en la tabla se aprecia un crecimiento paralelo entre las dos series, con un error relativamente bajo (del orden de 0.2 a 0.3 unidades). Sin embargo, se identifica un sesgo sistemático de sobreestimación: en todos los puntos el pronóstico se ubica ligeramente por encima del valor observado. Esto sugiere que, aunque el modelo es consistente en capturar la dinámica de las variables, tiende a generar proyecciones algo más optimistas que la realidad. El pronóstico para el cierre de 2024 y enero de 2025 mantiene esta tendencia. La proyección ubica el nivel esperado en torno a 145 unidades, mientras que el valor observado más reciente se aproxima a 144.9. La brecha es moderada y no desvirtúa el modelo, pero sí confirma la necesidad de tener cautela en la interpretación, dado que los escenarios futuros podrían estar ligeramente sobrevalorados. Desde el punto de vista empresarial y sectorial, los resultados son relevantes. Para las empresas cafeteras y actores del mercado, el modelo ofrece una herramienta útil para anticipar tendencias de precios y producción, facilitando la planificación de exportaciones, costos y estrategias de cobertura. No obstante, el sesgo de sobreestimación implica que confiar ciegamente en las cifras podría inducir decisiones demasiado optimistas en materia de inversión o manejo de inventarios. En consecuencia, se recomienda utilizar este pronóstico como guía de tendencia, complementándolo con análisis de sensibilidad y monitoreo de choques externos (climáticos, logísticos o de demanda internacional) que impactan de forma decisiva al sector cafetero colombiano.

Pronóstico del modelo automático SARIMA fuera de muestra, es decir, en enero 2025

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

Los resultados permiten anticipar que el sector cafetero colombiano enfrentará un inicio de 2025 con precios de exportación relativamente altos, lo que representa una ventaja competitiva para empresas exportadoras. Sin embargo, el sesgo identificado sugiere que el mercado real podría ubicarse en niveles levemente inferiores, por lo que es prudente interpretar el pronóstico como una guía de tendencia y no como un valor absoluto.

Recomendaciones Estratégicas Basadas en las señales extraídas (STL y correlaciones)

Tendencia positiva del PECAFE:

Aprovechar la fase de crecimiento para cerrar contratos forward o ventas anticipadas a precios altos.

Incrementar la planeación de exportaciones para el primer trimestre de 2025, asegurando ingresos.

Estacionalidad marcada en PNCAFE y PECAFE:

Planificar inventarios y flujos de caja en función de los meses de mayor volatilidad.

Coordinar la logística y almacenamiento para enfrentar picos de producción y precios.

IPC en crecimiento sostenido:

Reajustar estructuras de costos (insumos, salarios, transporte) para mitigar el efecto inflacionario.

Negociar contratos de suministro con cláusulas de indexación al IPC, protegiendo márgenes.

Ajuste de estrategias según el pronóstico de corto plazo

Si los precios se mantienen al alza (como proyecta el modelo):

Acelerar exportaciones y fortalecer posiciones en mercados internacionales.

Implementar políticas de cobertura financiera para asegurar rentabilidad en caso de correcciones de mercado.

Si el mercado real se ubica ligeramente por debajo de lo pronosticado (sesgo de sobreestimación):

Evitar decisiones de inversión excesivamente optimistas basadas solo en el pronóstico.

Usar análisis de sensibilidad y escenarios alternativos para gestionar riesgos.

Recomendación general:

Utilizar los pronósticos como herramientas de planeación estratégica, pero siempre complementarlos con monitoreo de factores externos (clima, logística, geopolítica, demanda internacional), que impactan de forma decisiva al sector cafetero.

# 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   144.2030
## 2 2024.833   144.4865
## 3 2024.917   145.0330
## 4 2025.000   146.0644
# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))
## [1] "Pronóstico para enero 2025: 2025 = 146.064423028787"