Series de tiempo: Extracción de señales y Modelo ARIMA
Daniel Medina
2025-03-12
knitr::opts_chunk$set(echo = TRUE, results = "show", message = FALSE, warning = FALSE)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",
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"))library(dplyr)
# Seleccionar variables
datos_cafe <- data_col %>%
select(PNCAFE, PECAFE, XCAF) %>%
na.omit()matriz_cor <- cor(datos_cafe)
matriz_cor## PNCAFE PECAFE XCAF
## PNCAFE 1.00000000 -0.05765602 0.2381612
## PECAFE -0.05765602 1.00000000 0.8403215
## XCAF 0.23816123 0.84032145 1.0000000# =========================
# 1. LIBRERÍAS
# =========================
library(dplyr)
library(corrplot)
# =========================
# 2. PREPARAR DATOS
# =========================
datos_cafe <- data_col %>%
select(PNCAFE, PECAFE, XCAF) %>%
na.omit()
# =========================
# 3. MATRIZ DE CORRELACIÓN
# =========================
matriz_cor <- cor(datos_cafe)
# Ver matriz numérica
print(matriz_cor)## PNCAFE PECAFE XCAF
## PNCAFE 1.00000000 -0.05765602 0.2381612
## PECAFE -0.05765602 1.00000000 0.8403215
## XCAF 0.23816123 0.84032145 1.0000000# =========================
# 4. CORRELOGRAMA (TIPO IMAGEN PRO)
# =========================
corrplot(
matriz_cor,
method = "circle", # círculos
type = "lower", # solo triángulo inferior
diag = FALSE, # sin diagonal
tl.col = "black", # etiquetas negras
tl.srt = 45, # rotación de nombres
addCoef.col = "black", # números dentro
number.cex = 0.8, # tamaño números
col = colorRampPalette(c("blue", "white", "darkblue"))(200)
)
La matriz de correlación evidencia relaciones diferenciadas entre las variables del sector cafetero. En particular, se observa una fuerte correlación positiva (0.84) entre el precio externo del café (PECAFE) y las exportaciones (XCAF), lo que indica que el comportamiento exportador del sector depende significativamente de las condiciones del mercado internacional.
Por otro lado, la relación entre la producción de café (PNCAFE) y las exportaciones (0.24) es positiva pero débil, lo que sugiere que el incremento en la producción no se traduce de manera proporcional en mayores exportaciones, posiblemente debido a factores como el aumento de demanda interna, inventarios de excedentes de los meses de alta producción para nivelar las exportaciones en los meses de baja producción o restricciones logísticas.
Finalmente, la correlación entre la producción y el precio externo (-0.06) es prácticamente nula, lo que indica que en el corto plazo la producción no responde directamente a los cambios en los precios internacionales, evidenciando posibles rezagos estructurales en la oferta del sector - (el sector agrícola colombiano tiene el cultivo del café como uno e los principales ejes).
En conjunto, estos resultados sugieren que el sector cafetero colombiano presenta una alta dependencia de variables externas, particularmente del precio internacional, más que de su propia capacidad productiva.
library(ggplot2)
library(dplyr)
# Datos
datos_cafe <- data_col %>%
select(PNCAFE, PECAFE, XCAF) %>%
na.omit()
# =========================
# 1. Producción vs Exportaciones
# =========================
ggplot(datos_cafe, aes(x = PNCAFE, y = XCAF)) +
geom_point(color = "steelblue", alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
theme_minimal() +
labs(
title = "Producción vs Exportaciones de Café",
x = "Producción de café (PNCAFE)",
y = "Exportaciones de café (XCAF)"
)
# =========================
# 2. Precio vs Exportaciones (CLAVE)
# =========================
ggplot(datos_cafe, aes(x = PECAFE, y = XCAF)) +
geom_point(color = "darkgreen", alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
theme_minimal() +
labs(
title = "Precio Internacional vs Exportaciones",
x = "Precio externo (PECAFE)",
y = "Exportaciones (XCAF)"
)
# =========================
# 3. Producción vs Precio
# =========================
ggplot(datos_cafe, aes(x = PNCAFE, y = PECAFE)) +
geom_point(color = "purple", alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE, color = "black") +
theme_minimal() +
labs(
title = "Producción vs Precio Internacional",
x = "Producción (PNCAFE)",
y = "Precio externo (PECAFE)"
)
PASO INDISPENSABLE: Declarar la (s) variable (s) como serie (s) temporal (es):
Variable 1
# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(data_col$PNCAFE, start = c(2012, 1), frequency = 12)Variable 2
# Convertir/declarar el ISE en serie de tiempo mensual
variable2_ts <- ts(data_col$PECAFE, start = c(2012, 1), frequency = 12)Variable 3
# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable3_ts <- ts(data_col$XCAF, start = c(2012, 1), frequency = 12)Extracción de señales
Gráfico inicial de la variable 1 en niveles -Original
library(ggplot2)
library(plotly)
# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable1_ts)
# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)),
y = variable1)) +
geom_line(color = "grey", linewidth = 0.4) + # Cambiado 'size' por 'linewidth'
geom_point(color = "black", size = 0.1) +
ggtitle("Variable 1: Serie original- producción de cafe") +
xlab("Tiempo") +
ylab("Unidad Variable 1") +
theme_minimal()
ggplotly(grafico_serie)Extracción señales variable 1- produccion de cáfe
# 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- producción de cafe",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Descomposición de la Producción de Café
Tendencia
La tendencia muestra un crecimiento sostenido desde 2012 hasta aproximadamente 2016, seguido de un comportamiento relativamente estable con ligeras fluctuaciones hasta 2020. Posteriormente, se evidencia una caída en la producción alrededor de 2021–2023, posiblemente asociada a factores climáticos, logísticos o económicos. Hacia el final del periodo se observa una recuperación parcial, aunque sin alcanzar los niveles máximos previos.
Estacionalidad
La componente estacional presenta un patrón altamente repetitivo a lo largo del tiempo, lo que indica la existencia de ciclos productivos bien definidos en el sector cafetero. Este comportamiento es consistente con los periodos de cosecha del café en Colombia, donde se registran picos y caídas regulares durante el año.
Componente irregular
El componente irregular muestra fluctuaciones aleatorias sin un patrón definido, aunque con algunos picos atípicos que podrían estar asociados a choques externos como condiciones climáticas adversas, variaciones en costos o disrupciones en la cadena de suministro.
En conjunto, la producción de café en Colombia presenta una fuerte estacionalidad y una tendencia que refleja tanto periodos de crecimiento como de ajuste, lo que evidencia que el sector está influenciado por factores estructurales y coyunturales. La presencia de variaciones irregulares sugiere además la vulnerabilidad del sector ante shocks externos.
Extracción señales variable 2- precio externo del cáfe
# 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- precio externo del cafe ",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Descomposición del Precio Externo del Café
La descomposición de la serie temporal del precio externo del café permite identificar patrones relevantes en su comportamiento, evidenciando la influencia de factores internacionales en el sector cafetero colombiano.
Tendencia
La tendencia del precio muestra una disminución inicial entre 2012 y aproximadamente 2015, seguida de un periodo relativamente estable hasta 2019. A partir de este punto, se observa un crecimiento significativo, especialmente desde 2020 en adelante, alcanzando niveles máximos hacia el final del periodo analizado. Este comportamiento puede estar asociado a cambios en la oferta global, condiciones climáticas internacionales y variaciones en la demanda mundial.
Estacionalidad
La componente estacional presenta un patrón repetitivo a lo largo del tiempo, aunque con menor intensidad en comparación con la producción. Esto indica que, si bien existen ciclos dentro del año, el precio internacional está menos determinado por factores estacionales y más por condiciones del mercado global.
Componente irregular
El componente irregular muestra fluctuaciones importantes, especialmente en los periodos más recientes, lo que refleja la alta volatilidad del precio internacional del café. Estos cambios pueden estar asociados a choques externos como crisis económicas, cambios en la oferta global o eventos climáticos.
Interpretación global
En conjunto, el precio externo del café presenta una dinámica dominada por factores internacionales, con una tendencia creciente en los últimos años y una volatilidad significativa. Esto evidencia que el sector cafetero colombiano está altamente expuesto a las condiciones del mercado global.
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 - exportaciones de cafe",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Interpretación – Descomposición de las Exportaciones de Café
Tendencia
La tendencia de las exportaciones muestra un comportamiento relativamente estable entre 2012 y 2019, seguido de un crecimiento significativo a partir de 2020. Este incremento se intensifica hacia el final del periodo, alcanzando los niveles más altos observados. Este comportamiento es consistente con el aumento del precio internacional del café, lo que sugiere una fuerte influencia de factores externos sobre el desempeño exportador.
Estacionalidad
La componente estacional presenta un patrón claro y repetitivo, con fluctuaciones marcadas a lo largo del año. Esto indica que las exportaciones están fuertemente influenciadas por los ciclos productivos del café, reflejando la sincronización entre la cosecha y la actividad exportadora.
Componente irregular
El componente irregular evidencia variaciones importantes en ciertos periodos, con la presencia de valores atípicos que podrían estar asociados a choques externos, como cambios en la demanda internacional, problemas logísticos o condiciones climáticas que afectan la oferta.
Interpretación global
En conjunto, las exportaciones de café presentan una dinámica caracterizada por una fuerte estacionalidad y una tendencia creciente en los últimos años, lo que refleja tanto la capacidad productiva del sector como la influencia de las condiciones del mercado internacional. La volatilidad observada en el componente irregular sugiere además una alta sensibilidad a factores externos.
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 Producción de café
# 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)Al comparar la serie original con la serie ajustada por estacionalidad, se observa que la variabilidad de la producción de café está fuertemente influenciada por patrones estacionales asociados a los ciclos de cosecha. Al eliminar este componente, la serie ajustada permite identificar con mayor claridad la tendencia subyacente, evidenciando un comportamiento más estable. Esto indica que los cambios en la producción no responden únicamente a factores económicos, sino principalmente a dinámicas propias del ciclo productivo agrícola.
Gráfico serie original VS ajustada Variable 2 - Precio externo del café colombiano
# 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)Al comparar la serie original con la ajustada por estacionalidad, se observa que ambas presentan un comportamiento muy similar, lo que indica que el precio externo del café no está fuertemente influenciado por patrones estacionales. En cambio, su dinámica está dominada por cambios en la tendencia y fluctuaciones asociadas a factores del mercado internacional, evidenciando una alta volatilidad y dependencia de condiciones externas.
Gráfico serie original VS ajustada Variable 3 exportaciones
# 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)Al comparar la serie original con la ajustada por estacionalidad, se observa que, aunque existen fluctuaciones periódicas, la dinámica de las exportaciones está principalmente determinada por la tendencia creciente en los últimos años. La serie ajustada permite evidenciar con mayor claridad este crecimiento, mostrando que el comportamiento exportador responde tanto a ciclos productivos como a factores externos, especialmente el aumento del precio internacional.
Ahora graficamos serie original vs tendencia
Primero se debe obtener la tendencia de cada variable y luego graficarla
Tendencia Variable 1 produccion de cafe
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("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)La tasa de crecimiento de la producción de café presenta una alta volatilidad en la serie original, evidenciando fluctuaciones significativas en el corto plazo. Sin embargo, al analizar la tendencia, se observa un crecimiento sostenido hasta aproximadamente 2016, seguido de un periodo de estabilidad y posterior desaceleración entre 2020 y 2022. Hacia el final del periodo se identifica una recuperación en el crecimiento, lo que sugiere un ajuste positivo reciente en la dinámica productiva del sector.
La diferencia entre la serie original y la tendencia indica que el crecimiento del sector está sujeto a choques de corto plazo, pero mantiene una dinámica estructural definida en el tiempo.
Tendencia Variable 2 precio externo del cafe
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("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)La serie del precio externo del café presenta fluctuaciones en el corto plazo; sin embargo, la tendencia permite identificar una dinámica clara: una caída inicial, seguida de un periodo de estabilidad y, posteriormente, un crecimiento sostenido a partir de 2020. Este comportamiento evidencia que el precio internacional está influenciado principalmente por factores globales y muestra una fase reciente de fuerte expansión.
Tendencia Variable 3 exportaciones
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)La serie de exportaciones presenta alta variabilidad en el corto plazo; sin embargo, la tendencia evidencia un comportamiento relativamente estable hasta 2020, seguido de un crecimiento marcado en los últimos años. Este aumento coincide con el incremento del precio internacional del café, lo que sugiere que el dinamismo reciente de las exportaciones está impulsado principalmente por factores externos.
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] 156print(length(tasa_crecimiento_var1))## [1] 156print(length(tasa_tendencia_var1))## [1] 156*Gráfico variable original y tendencia variable 1: tasa de crecimiento anual - produccion de cafe**
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)La tasa de crecimiento de la producción de café presenta una alta volatilidad en el corto plazo, con variaciones significativas entre periodos positivos y negativos. No obstante, la tendencia permite identificar un ciclo claro: un crecimiento elevado al inicio, seguido de una desaceleración progresiva hasta valores cercanos o negativos, y una recuperación reciente antes de una nueva caída al final del periodo. Esto sugiere que la dinámica productiva del sector es cíclica y altamente sensible a choques temporales.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2 - precio externo del cafe
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))
# Verificar longitudes
print(length(fechas_corregidas_var2))## [1] 156print(length(tasa_crecimiento_var2))## [1] 156print(length(tasa_tendencia_var2))## [1] 156# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable2: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)La tasa de crecimiento del precio externo del café presenta una alta volatilidad, con fuertes expansiones y contracciones a lo largo del tiempo. Sin embargo, la tendencia evidencia ciclos marcados: periodos de crecimiento acelerado seguidos de caídas pronunciadas y posteriores recuperaciones. Este comportamiento refleja la naturaleza cíclica y altamente sensible del precio internacional a factores del mercado global. A diferencia de la producción, el precio no sigue una dinámica estable, sino ciclos pronunciados, lo que explica su fuerte impacto sobre las exportaciones.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 3 - exportaciones de cáfe
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_var3 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var3))
# Verificar longitudes
print(length(fechas_corregidas_var3))## [1] 156print(length(tasa_crecimiento_var3))## [1] 156print(length(tasa_tendencia_var3))## [1] 156# Gráfico de la tasa de crecimiento anual variable 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)La tasa de crecimiento de las exportaciones de café presenta una alta volatilidad, con variaciones significativas entre periodos de expansión y contracción. No obstante, la tendencia permite identificar ciclos claros: fases de crecimiento moderado, caídas pronunciadas y una recuperación reciente con tasas positivas. Este comportamiento sugiere que el dinamismo exportador está influenciado por factores externos, especialmente el precio internacional, y no sigue una trayectoria de crecimiento estable.
La volatilidad en el crecimiento de las exportaciones confirma que el sector cafetero colombiano depende más de condiciones del mercado internacional que de una dinámica productiva estable.
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 3 que es la elegida para pronosticar
El código siguiente divide una serie temporal (variable3_ts) en dos subconjuntos:
Conjunto de entrenamiento (train): Datos desde enero de 2012 hasta septiembre de 2025. Conjunto de prueba (test): Datos desde octubre de 2025 hasta diciembre de 2025.
Esto se hace para evaluar el desempeño de modelos de predicción en datos no vistos.
# Esta división idealmente podria se 80%-70% de los datos para entrenamiento y 20%-30% para prueba o test
# En este ejemplo el conjunto de entrenamiento es: Enero 2012-Septiembre 2025 y el conjunto de prueba o test: noviembre 2025-diciembre 2025
train_size <- length(variable3_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable3_ts, end = c(2025, 9)) # Entrenamiento hasta septiembre 2025
test_ts <- window(variable3_ts, start = c(2025, 10)) # Prueba inicia desde oct2025Modelo 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(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.8485 0.3959 -0.5332
## s.e. 0.0628 0.0806 0.0683
##
## sigma^2 = 2.275e+09: log likelihood = -1998.25
## AIC=4004.5 AICc=4004.75 BIC=4016.9
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3932.893 47119.03 33776.97 -1.850663 14.93934 0.5269042 0.01120106Estimació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.848507 0.062776 -13.5164 < 2.2e-16 ***
## ma1 0.395889 0.080598 4.9119 9.02e-07 ***
## ma2 -0.533167 0.068283 -7.8082 5.80e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Ajuste del modelo ARIMA(1,1,2) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts,
order = c(1, 1, 2)) # Especificamos directamente (p=4, d=1, q=2)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.8485 0.3959 -0.5332
## s.e. 0.0628 0.0806 0.0683
##
## sigma^2 = 2.275e+09: log likelihood = -1998.25
## AIC=4004.5 AICc=4004.75 BIC=4016.9
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3932.893 47119.03 33776.97 -1.850663 14.93934 0.5269042 0.01120106Validació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(1,1,2)
## Q* = 34.827, df = 21, p-value = 0.0295
##
## Model df: 3. Total lags used: 24Pronó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 interactivoInterpretación modelo automatico (1,1,2):El modelo automático (4,1,2) 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 2025.750 547763.4 493522.0
## 2 2025.833 592944.4 501971.2
## 3 2025.917 524602.7 494802.0Ahora pronosticamos con el modelo automatico fuera del periodo de análisis, es decir enero 2026
Es decir, le sumamos al periodo de prueb auna observación más. Es decir, se estan pronosticando 4 observaciones o meses.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo trimestre
next_month_forecast_auto <- data.frame(
Tiempo = time(next_forecast_auto$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast_auto$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast_auto)## Tiempo Pronostico
## 1 2025.750 493522.0
## 2 2025.833 501971.2
## 3 2025.917 494802.0
## 4 2026.000 500885.1# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_auto, 1)
print(paste("Pronóstico para enero 2026:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2026: 2026 = 500885.115420756"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,0,1)[12]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## -0.8728 0.4146 -0.5324 0.1545
## s.e. 0.0544 0.0762 0.0676 0.0753
##
## sigma^2 = 2.23e+09: log likelihood = -1996.18
## AIC=4002.37 AICc=4002.74 BIC=4017.86A continuación, se crea el objeto darima para luegO poder graficar los valores reales y observados:
# Cargar el paquete necesario
library(forecast)
# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12] #Modelo identificado en el paso anterior
darima <- Arima(train_ts,
order = c(0, 1, 1), # (p,d,q) -> (0,1,1)
seasonal = list(order = c(1, 0, 0), # (P,D,Q) -> (1,0,0)
period = 12)) # Periodicidad estacional de 12 meses
# Mostrar resumen del modelo ajustado
summary(darima)## Series: train_ts
## ARIMA(0,1,1)(1,0,0)[12]
##
## Coefficients:
## ma1 sar1
## -0.5197 0.1971
## s.e. 0.0705 0.0805
##
## sigma^2 = 2.262e+09: log likelihood = -1998.35
## AIC=4002.69 AICc=4002.84 BIC=4011.99
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 3615.563 47129.43 35382.69 -1.910218 15.66074 0.5519526 0.02438643Validación de residuales del modelo automatico SARIMA
En el correlograma de residuos siguiente se observa que, mejora la correlación de los residuos frente a lso dos modelos anteriores. Sin embargo, al comparar los valores reales VS pronosticados se determina una poca coincidencia. Sigue funcionando mejor el modelo automatico (4,1,2)
# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(1,0,0)[12]
## Q* = 25.302, df = 22, p-value = 0.2829
##
## Model df: 2. Total lags used: 24Pronó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 interactivoPronóstico del modelo automático SARIMA en el set de prueba-Tabla
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast <- forecast(auto_arima_model, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
Tiempo = time(arima_forecast$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table)## Tiempo Observado Pronosticado
## 1 2025.750 547763.4 502617.4
## 2 2025.833 592944.4 506938.1
## 3 2025.917 524602.7 520137.1Pronó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.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
Tiempo = time(next_forecast$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast)## Tiempo Pronostico
## 1 2025.750 502617.4
## 2 2025.833 506938.1
## 3 2025.917 520137.1
## 4 2026.000 529353.7# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2025: 2026 = 529353.664622314"Conclusión:
El modelo automático ARIMA(4,1,2) fue el que mejor desempeño mostró en la comparación entre los datos reales y los pronosticados dentro del periodo de prueba (oct.nov.dic2024). Destacó por su mayor precisión en la captura de los puntos de quiebre, lo que lo hace el más confiable.
No obstante, al analizar los residuos de los modelos, se identifican posibles áreas de mejora para robustecer los pronósticos en los tres casos evaluados. Algunas estrategias podrían incluir la aplicación de una transformación logarítmica o trabajar desde el inicio con la serie ajustada por estacionalidad.
By Julieth Cerón