Series de tiempo: Extracción de señales y Modelo ARIMA para la producción del café
Pablo López & Natalya López
2025-03-23
Instalar/Cargar librerias necesarias para el análisis
#Cargar librerías necesarias
library(readxl) # Para leer archivos Excel
library(tseries) # Para pruebas de estacionariedad
library(forecast) # Para modelado ARIMA y pronósticos
library(ggplot2) # Para visualización de datos
library(plotly) # Para gráficos interactivos
library(timetk) Cargar base de datos
library(readxl)
library(readxl)
Cafe <- read_excel("Cafe.xlsx", col_types = c("date",
"numeric", "numeric", "numeric"))
View(Cafe)PASO INDISPENSABLE: Declarar la (s) variable (s) como serie (s) temporal (es):
Variable 1: Producción del café
# Convertir/declarar variable 1=ENER en serie de tiempo mensual
pncafe_ts <- ts(Cafe$PNCAFE, start = c(2012, 1), frequency = 12)
pncafe_ts## Jan Feb Mar Apr May Jun Jul
## 2012 535.0000 571.0000 576.0000 580.0000 689.0000 714.0000 668.0000
## 2013 877.0000 625.0000 617.0000 970.0000 937.0000 913.0000 1031.0000
## 2014 1011.0000 874.0000 828.0000 832.0000 1050.0000 944.0000 1236.0000
## 2015 1088.0000 1029.0000 800.0000 924.0000 1165.0000 1240.0000 1463.0000
## 2016 1136.0000 1096.0000 944.0000 1043.0000 1163.0000 1158.0000 1102.0000
## 2017 1275.0000 1293.0000 1020.0000 834.0000 901.0000 1049.0000 1373.0000
## 2018 1131.0000 1212.0000 1037.0000 874.0000 1188.0000 1087.0000 1051.0000
## 2019 1296.0000 1106.0000 914.0000 1031.0000 1115.0000 1211.0000 1317.0000
## 2020 1050.0000 1001.0000 806.0000 744.0000 1186.0000 1362.0000 1310.0000
## 2021 1081.0000 1107.0000 1050.0000 810.0000 616.0000 1052.0000 1209.0000
## 2022 868.0000 928.0000 914.0000 750.0000 1017.0000 951.0000 943.5992
## 2023 868.0000 1025.0000 799.0000 565.8672 806.1955 956.0000 947.0000
## 2024 959.0000 961.0000 865.5326 742.0000 1120.0000 1172.0000 1158.8838
## Aug Sep Oct Nov Dec
## 2012 565.0000 519.0000 653.0000 770.0000 904.0000
## 2013 770.0000 860.0000 1058.0000 1113.0000 1115.0000
## 2014 1151.0000 912.0000 1101.0000 1115.0000 1086.0000
## 2015 1264.0000 1058.0000 1368.0000 1322.0000 1454.0000
## 2016 1189.0000 1034.0000 1395.0000 1653.0000 1319.0000
## 2017 1294.0000 1228.0000 1073.0000 1304.0000 1550.0000
## 2018 1258.0000 1050.0000 1086.0000 1300.0000 1283.0000
## 2019 1119.0000 1088.0000 1369.0000 1506.0000 1680.0000
## 2020 1091.0000 995.0000 1159.0000 1443.0000 1743.0000
## 2021 915.0000 1209.0000 1012.0000 1131.0000 1385.0000
## 2022 949.0000 834.0000 888.0000 1060.0000 981.0000
## 2023 872.0000 849.0000 1157.4585 1282.1066 1220.0000
## 2024 1049.0000 1071.2334 1339.1550 1761.4095 1798.2306Variable 2: Precio externo del café colombiano
# Convertir/declarar el ISE en serie de tiempo mensual
PECAFE_ts <- ts(Cafe$PECAFE, start = c(2012, 1), frequency = 12)
PECAFE_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2012 256.1168 246.3010 226.0661 215.2930 209.5097 185.8843 200.9219 188.5771
## 2013 168.8300 162.5343 161.7474 162.1890 160.0019 148.4153 147.5277 144.2926
## 2014 131.9742 159.6775 210.1481 216.0827 216.4171 195.8677 193.0277 210.4074
## 2015 185.3126 177.0854 155.8761 156.4597 152.3323 151.1303 145.9206 147.1319
## 2016 136.4284 137.1914 144.3697 144.5490 144.5803 153.9430 164.0368 160.7271
## 2017 163.7090 164.8057 159.2329 156.0010 151.2374 146.8607 150.9077 156.4565
## 2018 143.4810 142.3511 139.7903 138.9967 140.5035 139.1723 134.8471 130.9932
## 2019 128.6758 128.9696 125.4671 124.0727 123.1348 132.7700 138.4145 129.8942
## 2020 150.7962 144.1289 157.3050 164.6219 155.8180 149.2377 151.2900 167.6310
## 2021 175.2395 177.2779 179.1478 181.6810 196.2590 206.1064 214.4400 226.9009
## 2022 293.1955 305.5184 287.9070 292.5780 285.1871 301.2423 286.9865 292.5200
## 2023 218.2400 237.0400 228.2448 232.3695 229.0759 215.4367 190.3325 188.7070
## 2024 206.5600 210.1300 207.5470 237.3114 232.3527 249.9400 257.9800 261.0000
## Sep Oct Nov Dec
## 2012 189.2787 183.1594 171.1003 166.0471
## 2013 139.1247 135.1697 125.3610 125.9203
## 2014 208.0950 223.1000 206.7437 193.3329
## 2015 136.4220 142.7784 138.4003 139.6919
## 2016 167.6153 170.5529 179.7767 160.9187
## 2017 151.5747 145.0094 143.8483 142.5123
## 2018 126.5253 138.8045 140.7503 129.6229
## 2019 131.3550 131.8783 143.5175 160.1652
## 2020 169.8400 156.7536 161.6000 169.9900
## 2021 238.3200 257.1257 273.6352 291.9477
## 2022 296.4571 269.4900 225.1000 223.8900
## 2023 186.0000 184.9700 194.4567 206.9235
## 2024 276.2520 280.0396 293.9385 340.5167Variable 3: Exportaciones de café
# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
XCAF_ts <- ts(Cafe$XCAF, start = c(2012, 1), frequency = 12)
XCAF_ts## Jan Feb Mar Apr May Jun Jul Aug
## 2012 197295.8 186693.3 203636.0 121442.5 167643.9 162433.1 151832.1 180244.0
## 2013 177150.5 171569.2 127181.2 152040.1 168982.2 138026.6 159759.0 169039.9
## 2014 142633.8 202636.1 153516.9 224406.6 213123.9 161247.5 229475.0 218584.5
## 2015 288671.5 251209.4 211129.0 182545.4 194717.5 160584.9 239676.4 218630.5
## 2016 163588.4 184225.8 179113.8 281209.6 152586.3 176288.8 128997.8 169424.1
## 2017 232188.4 214883.5 295930.9 217795.8 196637.3 137977.0 230431.8 211174.9
## 2018 255058.1 187220.6 201578.7 205465.6 190580.9 147559.5 166218.3 213792.4
## 2019 260699.1 221412.3 193807.5 186052.9 143552.4 147796.7 205010.7 218999.2
## 2020 249836.3 215605.6 170321.5 157384.0 141864.0 194805.8 270277.4 201786.8
## 2021 255455.2 267999.8 278224.9 243363.4 124318.7 115655.4 309123.1 322175.0
## 2022 332201.4 383477.5 412235.8 327940.2 275625.7 369007.8 394622.5 323525.4
## 2023 250131.8 261153.3 292938.5 212748.0 251607.5 225461.4 227622.0 226253.7
## 2024 231800.6 258106.5 238506.5 261112.5 268005.6 285915.9 291114.4 318115.4
## Sep Oct Nov Dec
## 2012 123549.6 159793.2 144269.2 166159.5
## 2013 152688.2 152603.0 173174.7 191779.4
## 2014 216896.0 252496.3 267402.2 244052.0
## 2015 250151.3 179340.6 191276.3 217993.9
## 2016 196074.9 193648.9 219217.4 429058.8
## 2017 240608.9 215414.3 198250.7 203778.0
## 2018 179237.1 194640.1 191052.0 216342.7
## 2019 167456.2 195326.7 193952.8 242374.7
## 2020 217897.4 173371.9 227635.8 316675.7
## 2021 268316.6 327953.2 316349.7 376919.2
## 2022 339103.7 323641.1 293141.5 348771.0
## 2023 234930.5 188621.9 251091.2 306637.2
## 2024 284093.6 338376.9 325234.9 462940.1Extracción de señales
*Muchas series de tiempo son una combinación de varias influencias. Es por eso que, separar la tendencia, la estacionalidad y los componentes aleatorios permite entender mejor qué está impulsando los cambios en la serie.
Si analizamos la producción del café, podríamos querer saber si el crecimiento se debe a una tendencia real o a fluctuaciones estacionales.
- Los modelos de pronóstico funcionan mejor cuando las señales subyacentes están bien definidas. Por ejemplo, si eliminamos la estacionalidad de la producción del café, los modelos predictivos pueden enfocarse en la tendencia real y reducir errores.
*Detectar cambios inesperados en la serie es más fácil cuando se eliminan componentes predecibles. Ejemplo: Si hay una caída abrupta en la variable podemos verificar si es una anomalía (ruido) o un cambio estructural en la producción del café.
En conclusión, la descomposición de series de tiempo permite comprender mejor los datos, mejorar predicciones y tomar decisiones más estratégicas. Es una herramienta clave en la analítica de negocios, especialmente en entornos donde las fluctuaciones en los datos pueden afectar inversiones, políticas económicas y estrategias empresariales.
Gráfico inicial de la variable 1: Producción del café en niveles -Original
library(ggplot2)
library(plotly)
# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
Cafe$PNCAFE <- as.numeric(pncafe_ts)
# Crear el gráfico
grafico_serie <- ggplot(Cafe, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(Cafe)),
y = PNCAFE)) +
geom_line(color = "grey", linewidth = 0.4) + # Cambiado 'size' por 'linewidth'
geom_point(color = "black", size = 0.1) +
ggtitle("Producción del café: Serie original") +
xlab("Tiempo") +
ylab("Toneladas Producción del café") +
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_pncafe <- stl(pncafe_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_pncafe <- data.frame(
Time = rep(time(pncafe_ts), 4), # Tiempo repetido para cada componente (son 4 componentes)
Value = c(stl_decomp_pncafe$time.series[, "seasonal"],
stl_decomp_pncafe$time.series[, "trend"],
stl_decomp_pncafe$time.series[, "remainder"],
pncafe_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(pncafe_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_pncafe, 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 Producción del café",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Interpretación Descomposición temporal para producción del café
Componente estacional
Las caidas o incrementos de la producción del café estacional ocurren a intervalos similares (Primer trimestre y cuarto trimestre de cada año), lo que indica que el comportamiento estacional se mantiene a lo largo del tiempo.
Así, el componente estacional sugiere que hay fluctuaciones repetitivas en la producción del café a lo largo de cada año. Comprender estos patrones permite ajustar estrategias de apoyo a la metodologia de producción y realizar pronósticos más precisos para futuras decisiones en aquellos periodos que se ve afectado.
Extracción señales variable 2
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp_pecafe <- stl(PECAFE_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_pecafe <- data.frame(
Time = rep(time(PECAFE_ts), 4), # Tiempo repetido para cada componente
Value = c(stl_decomp_pecafe$time.series[, "seasonal"],
stl_decomp_pecafe$time.series[, "trend"],
stl_decomp_pecafe$time.series[, "remainder"],
PECAFE_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(PECAFE_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_pecafe, aes(x = Time, y = Value, color = Component)) +
geom_line() +
facet_wrap(~Component, scales = "free_y", ncol = 1) +
theme_minimal() +
labs(title = "Descomposición temporal del Precio externo del café colombiano",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Extracción señales variable 3
# Cargar librerías necesarias
library(ggplot2)
library(plotly)
# Descomposición de la serie temporal
stl_decomp_XCAF <- stl(XCAF_ts, s.window = "periodic")
# Convertir la descomposición a un data frame para graficar con ggplot2
stl_df_xcaf <- data.frame(
Time = rep(time(XCAF_ts), 4), # Tiempo repetido para cada componente
Value = c(stl_decomp_XCAF$time.series[, "seasonal"],
stl_decomp_XCAF$time.series[, "trend"],
stl_decomp_XCAF$time.series[, "remainder"],
XCAF_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each = length(XCAF_ts))
)
# Crear gráfico con ggplot2
p <- ggplot(stl_df_xcaf, 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 las Exportaciones de café",
x = "Tiempo",
y = "Valor")
# Convertir a gráfico interactivo con plotly
ggplotly(p)Después de la descomposición temporal de cada variable, se extrae la variable ajustada por estacionalidad para graficarla junto con la serie original:
Se crea la variable1 ajustada por estacionalidad
# Extraer los componentes de la descomposición
pncafe_sa <- pncafe_ts - stl_decomp_pncafe$time.series[, "seasonal"]Se crea la variable2 ajustada por estacionalidad
# Extraer los componentes de la descomposición
pecafe_sa <- PECAFE_ts - stl_decomp_pecafe$time.series[, "seasonal"]Se crea la variable3 ajustada por estacionalidad
# Extraer los componentes de la descomposición
XCAF_sa <- XCAF_ts - stl_decomp_XCAF$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_pncafe <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(pncafe_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_pncafe <- ggplot() +
geom_line(aes(x = fechas_pncafe, y = pncafe_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_pncafe, y = pncafe_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Producción del café:Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Toneladas de Producción de café") +
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_pncafe)En esta ocasión se observa que en la serie ajustada por estacionalidad se pueden apreciar con mayor precisión cambios estructurales, facilitando el análisis sin la interferencia de variaciones temporales recurrentes.
Estos análisis permiten diseñar estrategias de apoyo al emprendimiento basadas en datos reales, sin sesgos estacionales. Además, proporciona información clave para la toma de decisiones estratégicas en el sector de la producción del café y, mejora la precisión de las estimaciones futuras al eliminar patrones estacionales que podrían distorsionar los análisis.
En conclusión, la descomposición de la serie temporal no solo facilita la identificación de tendencias de largo plazo, sino que también permite tomar decisiones más informadas, anticipando oportunidades y riesgos en la producción del café
Gráfico serie original VS ajustada Variable 2
# Crear vector de fechas correctamente alineado con la serie
fechas_pecafe <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(PECAFE_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_pecafe <- ggplot() +
geom_line(aes(x = fechas_pecafe, y = PECAFE_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_pecafe, y = pecafe_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Precio externo del café colombiano:Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Precio externo del café colombiano") +
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_pecafe)Gráfico serie original VS ajustada Variable 3
# Crear vector de fechas correctamente alineado con la serie
fechas_XCAF <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(XCAF_ts))
# Gráfico mejorado con fechas en el eje X
grafico_ajustada_XCAF <- ggplot() +
geom_line(aes(x = fechas_XCAF, y = XCAF_ts), color = "grey", size = 0.5, linetype = "solid", name = "Serie Original") +
geom_line(aes(x = fechas_XCAF, y = XCAF_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
ggtitle("Exportaciones de café:Serie Original vs Serie Ajustada por Estacionalidad") +
xlab("Tiempo") +
ylab("Exportaciones de café") +
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_XCAF)Ahora graficamos serie original vs tendencia
La extracción de la tendencia permite centrarse en los cambios estructurales de la serie.
Analizar la tendencia ayuda a prever escenarios futuros y anticipar posibles crisis o oportunidades en el sector o variable de análisis
Primero se debe obtener la tendencia de cada variable y luego graficarla
Tendencia Variable 1
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
pncafe_vec <- as.numeric(pncafe_ts)
tendencia_pncafe <- as.numeric(stl_decomp_pncafe$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(pncafe_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_pncafe <- ggplot() +
geom_line(aes(x = fechas, y = pncafe_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_pncafe, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Producción del café: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Toneladas de Producción del café") +
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_pncafe)Tendencia Variable 2
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
pecafe_vec <- as.numeric(PECAFE_ts)
tendencia_pecafe <- as.numeric(stl_decomp_pecafe$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(PECAFE_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_pecafe <- ggplot() +
geom_line(aes(x = fechas, y = pecafe_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_pecafe, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Precio externo del café colombiano: 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_pecafe)Tendencia Variable 3
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
xcaf_vec <- as.numeric(XCAF_ts)
tendencia_xcaf <- as.numeric(stl_decomp_XCAF$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(XCAF_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_xcaf <- ggplot() +
geom_line(aes(x = fechas, y = xcaf_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_xcaf, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Exportaciones de café: 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_xcaf)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_pncafe <- (pncafe_ts[(13:length(pncafe_ts))] / pncafe_ts[1:(length(pncafe_ts) - 12)] - 1) * 100
tasa_tendencia_pncafe <- (tendencia_pncafe[(13:length(tendencia_pncafe))] / tendencia_pncafe[1:(length(tendencia_pncafe) - 12)] - 1) * 100
# Crear vector de fechas corregido, es decir que inicie desde enero 2013
fechas_corregidas_pncafe <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_pncafe))
# Verificar longitudes
print(length(fechas_corregidas_pncafe))## [1] 144print(length(tasa_crecimiento_pncafe))## [1] 144print(length(tasa_tendencia_pncafe))## [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_pncafe <- ggplot() +
geom_line(aes(x = fechas_corregidas_pncafe, y = tasa_crecimiento_pncafe), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_pncafe, y = tasa_tendencia_pncafe), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Producción de café: 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_pncafe)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_pecafe <- (PECAFE_ts[(13:length(PECAFE_ts))] / PECAFE_ts[1:(length(PECAFE_ts) - 12)] - 1) * 100
tasa_tendencia_pecafe <- (tendencia_pecafe[(13:length(tendencia_pecafe))] / tendencia_pecafe[1:(length(tendencia_pecafe) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_pecafe <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_pecafe))
# Verificar longitudes
print(length(fechas_corregidas_pecafe))## [1] 144print(length(tasa_crecimiento_pecafe))## [1] 144print(length(tasa_tendencia_pecafe))## [1] 144# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_pecafe <- ggplot() +
geom_line(aes(x = fechas_corregidas_pecafe, y = tasa_crecimiento_pecafe), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_pecafe, y = tasa_tendencia_pecafe), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Precio externo del café colombiano: 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_pecafe)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_xcaf <- (XCAF_ts[(13:length(XCAF_ts))] / XCAF_ts[1:(length(XCAF_ts) - 12)] - 1) * 100
tasa_tendencia_xcaf <- (tendencia_xcaf[(13:length(tendencia_xcaf))] / tendencia_xcaf[1:(length(tendencia_xcaf) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_xcaf <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_xcaf))
# Verificar longitudes
print(length(fechas_corregidas_xcaf))## [1] 144print(length(tasa_crecimiento_xcaf))## [1] 144print(length(tasa_tendencia_xcaf))## [1] 144# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_xcaf <- ggplot() +
geom_line(aes(x = fechas_corregidas_xcaf, y = tasa_crecimiento_xcaf), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_xcaf, y = tasa_tendencia_xcaf), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Exportaciones de café: 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_xcaf)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
Un modelo ARIMA (Autoregressive Integrated Moving Average) es una herramienta estadística utilizada para analizar y predecir series de tiempo. En términos simples, es como una “bola de cristal matemática” que usa datos pasados para estimar valores futuros, especialmente útil en finanzas para preveer precios, ventas o ingresos.
Desglose del nombre ARIMA
- Autoregressive o autorregresivo (AR) → Usa valores pasados para predecir el futuro
- Integrated (I) → Ajusta tendencias en los datos para hacerlos estacionarios (sin patrones cambiantes en el tiempo).
3.MAving Average o media móvil (MA)- → Suaviza fluctuaciones aleatorias a partir de errores pasados.
ARIMA combina estos tres elementos para crear una predicción más precisa.
Metodología Box-Jenkins para aplicar un modelo ARIMA
La metodología Box-Jenkins es un enfoque sistemático para construir modelos ARIMA con el objetivo de analizar y pronosticar series de tiempo. Fue desarrollada por George Box y Gwilym Jenkins y se basa en cuatro etapas clave:
1️⃣ Identificación 2️⃣ Estimación 3️⃣ Validación 4️⃣ Pronóstico
Se usa especialmente cuando se quiere encontrar el modelo ARIMA más adecuado para una serie de tiempo.
Antes de empezar a aplicar la metodología BOX-JENKINS, lo ideal es dividir el conjunto de datos de prueba y entrenamiento
✅ Para entrenar el modelo con datos históricos sin usar información futura. ✅ Para evaluar la precisión del modelo comparando sus predicciones con los datos reales de prueba.
💡 Esta división es clave en modelos predictivos para evitar sobreajuste y evaluar el rendimiento en datos no vistos.
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.
pncafe_ts <- ((pncafe_ts))# 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(pncafe_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(pncafe_ts, end = c(2024, 9)) # Entrenamiento hasta septiembre 2024
test_ts <- window(pncafe_ts, start = c(2024, 10)) # Prueba inicia desde oct2024Paso 1: Identificación del modelo
Identificar estacionariedad
En el análisis de series de tiempo, una serie es estacionaria si su comportamiento es constante a lo largo del tiempo, es decir:
✅ Su promedio no cambia con el tiempo. ✅ Su variabilidad (qué tanto fluctúa) se mantiene estable. ✅ Su relación con valores pasados es siempre la misma.
¿Cómo saber si una serie es estacionaria?
1️⃣ Observando un gráfico Si el gráfico de la serie muestra una tendencia creciente o decreciente, o si las variaciones se hacen más grandes con el tiempo, la serie probablemente no es estacionaria.
2️⃣ Usando la prueba de Dickey-Fuller Aumentada (ADF) Es una prueba estadística que nos dice si la serie tiene una tendencia fuerte. Si el p-valor de la prueba es mayor a 0.05, significa que la serie no es estacionaria.
¿Cómo hacer una serie estacionaria?
Si encontramos que la serie no es estacionaria, podemos transformarla para que lo sea:
✅ Diferenciación: Restamos cada valor con su valor anterior. Esto elimina tendencias crecientes o decrecientes. ✅ Tomar el logaritmo: Si la variabilidad crece con el tiempo, aplicar un logaritmo estabiliza la varianza. ✅ Eliminar tendencias o ajustar estacionalidad: Si hay patrones repetitivos, podemos restarlos o modelarlos por separado.
Test de Dickey-Fuller
El test de Dickey-Fuller aumentado (ADF) se usa para verificar si una serie temporal es estacionaria, es decir, si sus propiedades estadísticas (media y varianza) permanecen constantes en el tiempo.
HO: Serie no estacionaria HI: Serie estacionaria
¿Qué significa el p-valor?
Si el p-valor es bajo (< 0.05) → Rechazamos la hipótesis nula y concluimos que la serie es estacionaria. Si el p-valor es alto (> 0.05) → No podemos rechazar la hipótesis nula, lo que indica que la serie no es estacionaria.
A continuación se aplica el test ADF para validar estacionariedad en el conjunto de entrenamiento de la variable 1, que es la elegida para pronosticar:
library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la Producción del café (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test##
## Augmented Dickey-Fuller Test
##
## data: train_ts
## Dickey-Fuller = -3.1273, Lag order = 5, p-value = 0.1066
## alternative hypothesis: stationaryEl test ADF en la variable 1 arrojó un p-value igual a 0.13, este valor es mayor a 0.05, por tanto la serie es no estacionaria. De ese modo se debe ejecutar el código siguiente para diferenciar una vez la variable 1 y luego volver a aplicar el test ADF a esa serie diferenciada una vez:
#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable 1 , una sola vez:
train_diff <- diff(train_ts, differences = 1) Diferenciación en niveles variable 1
A continuación, se realiza el gráfico de la serie original y diferenciada (una vez) de la variable 1 para ver graficamente el cambio o ajuste:
library(ggplot2)
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), pncafe = as.numeric(train_ts)), aes(x = Tiempo, y = pncafe)) +
geom_line(color = "blue") +
ggtitle("Producción del café:Serie Original") +
xlab("Tiempo") + ylab("Gwh")
library(ggplot2)
ggplotly(p2) # Convertir en gráfico interactivo # Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], pncafe_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = pncafe_Diff)) +
geom_line(color = "red") +
ggtitle("Producción de café:Serie Estacionaria (Una diferenciación)") +
xlab("Tiempo") + ylab("Producción de café diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoEjemplo Diferenciación en logaritmo
Cuando una serie de tiempo tiene una creciente varianza (lo que significa que la amplitud de las fluctuaciones aumenta con el tiempo), aplicar un logaritmo puede ayudar a estabilizar esta varianza. Muchas series económicas o financieras, como el precio de acciones o el Producto Interno Bruto (PIB), tienden a mostrar crecimiento exponencial o crecimiento en porcentaje (por ejemplo, tasas de crecimiento de doble dígito).
En conclusión, la aplicación de logaritmos en series de tiempo se realiza principalmente para lograr que la serie sea más estable, lineal y estacionaria. Esta transformación es relevante porque permite modelar mejor las series que siguen un crecimiento exponencial y facilita la aplicación de técnicas estadísticas que requieren estacionariedad.
A continuación se aplica la diferenciación logarítimica y la varible u objeto ahora se llama train_diff_log:
# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación
train_diff_log <- diff(log(train_ts), differences = 1)Ahora graficamos la serie orignal versus la serie diferenciada una vez con logaritmo
# Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), pncafe = as.numeric(train_ts)), aes(x = Tiempo, y = pncafe)) +
geom_line(color = "blue") +
ggtitle("pncafe:Serie Original") +
xlab("Tiempo") + ylab("pncafe")
ggplotly(p2) # Convertir en gráfico interactivo # Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], micro_Diff = as.numeric(train_diff_log)), aes(x = Tiempo, y = micro_Diff)) +
geom_line(color = "red") +
ggtitle("Producción de café:Serie Estacionaria (Una diferenciación en logaritmo)") +
xlab("Tiempo") + ylab("Producción de café diferenciada")
ggplotly(p3) # Convertir en gráfico interactivoAhora probamos estacionariedad en la serie diferenciada ( en nivel y logaritmo)
# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
adf_test_diff <- adf.test(train_diff)
print(adf_test_diff)##
## Augmented Dickey-Fuller Test
##
## data: train_diff
## Dickey-Fuller = -7.3148, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary# Segunda prueba de estacionariedad sobre la serie diferenciada en logaritmo
adf_test_diff_log <- adf.test(train_diff_log)
print(adf_test_diff_log)##
## Augmented Dickey-Fuller Test
##
## data: train_diff_log
## Dickey-Fuller = -7.3287, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryEn el test ADF se muestra que el valor que puede tomar d=1:
El p-value ya es menor a 0.05 con una primera diferencia en ambos casos: niveles o con logaritmo natural. Por tanto el valor que puede tomar d es igual a 1
Identificación manual de p y q
¿Qué hacen estos gráficos?
ACF (Autocorrelation Function)
Muestra la correlación de la serie con sus rezagos.
Ayuda a determinar el parámetro q en un modelo ARIMA(p, d, q).
PACF (Partial Autocorrelation Function)
Muestra la correlación parcial entre la serie y un rezago específico, eliminando el efecto de rezagos intermedios.
Ayuda a determinar el parámetro p en un modelo ARIMA(p, d, q).
Recordemos El eje X representa los rezagos (lags). El eje Y muestra la autocorrelación parcial en cada rezago. Las líneas azules punteadas son los intervalos de confianza (aproximadamente 95%). Si una barra sobrepasa estos límites, indica una autocorrelación significativa. Si las barras caen dentro de los límites, no son significativamente diferentes de cero
En el código siguiente se crean los correlogramas para determinar los posibles valores que puedeo tomar el parámetro p** y q:**
library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")
ggplotly(acf_plot)ggplotly(pacf_plot)Interpretación correlogramas
Se puede observar que los valores que podrian tomar p y q serian:
p=2* p=4 ; q=2* q=4 q=6 ; (P optimo=2) (q óptimo=2)
El modelo óptimo para esta variable seria (2,1,2)
Paso 2. Estimación manual del modelo
AIC y BIC: Se usan para comparar modelos; cuanto más bajos, mejor.
Métricas de evaluación
Mean Absolute Error (MAE)
Representa el error absoluto promedio entre las predicciones del modelo y los valores reales.
Root Mean Squared Error (RMSE)
Similar al MAE, pero da más peso a los errores grandes, porque eleva las diferencias al cuadrado antes de promediarlas.
Comparación con MAE: Como el RMSE es mayor que el MAE, es posible que haya algunos errores grandes que estén influyendo más en el RMSE.
Mean Absolute Percentage Error (MAPE) = 2.91%
Expresa el error en términos relativos, como porcentaje del valor real.
Interpretación: En promedio, el modelo se equivoca en un 2,91% al predecir la energía
Regla general: MAPE < 10% → Muy buen modelo ✅ 10%-20% → Modelo aceptable 👍 20%-50% → Modelo pobre ⚠️ 50% → Modelo muy malo ❌
En este caso, un MAPE de 2.91% sugiere un muy buen modelo para pronostico.
Estimación del modelo identificado (2,1,4)
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(2,1,4)) #Se va a estimar un modelo Arima de orden (2,1,4)
summary(manual_arima_model)## Series: train_ts
## ARIMA(2,1,4)
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4
## -0.5123 -0.6962 0.2398 0.2108 -0.4348 -0.5585
## s.e. 0.1542 0.0922 0.1358 0.0839 0.0714 0.0655
##
## sigma^2 = 26861: log likelihood = -988.64
## AIC=1991.28 AICc=1992.06 BIC=2012.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 14.32523 160.0987 124.9687 -0.3288151 12.23095 0.7462909
## ACF1
## Training set -0.02289033Significancia de coefientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.512290 0.154180 -3.3227 0.0008916 ***
## ar2 -0.696154 0.092213 -7.5494 4.373e-14 ***
## ma1 0.239824 0.135809 1.7659 0.0774146 .
## ma2 0.210780 0.083911 2.5119 0.0120068 *
## ma3 -0.434818 0.071404 -6.0896 1.132e-09 ***
## ma4 -0.558467 0.065496 -8.5267 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación significancia coeficientes
ar1, es altamente significativos (***), lo que significa que este coeficiente tiene un impacto importante en el modelo. Como al menos uno de los dos componentes e significativo, entonces continuo con la validación de residuos del modelo.
Paso 3. Validación de residuos del modelo estimado manual
La validación de residuos es crucial para determinar si el modelo ARIMA es adecuado o si necesita mejoras. El objetivo es verificar que los residuos (errores de predicción) se comporten como ruido blanco, es decir, sin patrones detectables.
1. Serie de residuos (gráfico superior)
Muestra cómo se comportan los errores a lo largo del tiempo. Idealmente, deberían oscilar alrededor de cero sin tendencias evidentes ni grandes acumulaciones de error. Problema posible: Se observan picos alrededor de 2020, lo que podría indicar cambios estructurales en los datos (como efectos de la pandemia en la demanda de energía que es la variable 1).
Función de Autocorrelación (gráfico inferior izquierdo, ACF de residuos)
Si el modelo es adecuado, los residuos no deben mostrar correlaciones significativas en el tiempo.
Interpretación: La mayoría de las barras están dentro de las líneas azules (intervalos de confianza). Sin embargo, algunos rezagos parecen salir del rango, lo que sugiere que aún puede haber estructura no capturada en los datos. Esto indica que el modelo podría mejorarse pero es un modelo aceptable.
Histograma de residuos con ajuste normal (gráfico inferior derecho)
Sirve para verificar si los errores siguen una distribución normal, lo cual es un supuesto clave en ARIMA. Interpretación: La curva roja representa la distribución normal teórica. Los residuos se acercan a la normalidad, pero hay algunos valores extremos (colas más gruesas de lo esperado). Esto indica que puede haber eventos atípicos o datos no bien explicados por el modelo (en este caso pandemia).
Conclusión y acciones recomendadas
✅ El modelo ARIMA(1,1,1)parece razonablemente bueno, pero tiene algunas señales de que podría mejorarse.
⚠️ Posibles mejoras:
Revisar la estructura del modelo: Se pueden probar otros órdenes ARIMA o incluso modelos más complejos como SARIMA o modelos con variables exógenas (ARIMAX).
Manejo de valores atípicos: Considerar incluir un término de intervención si eventos como la pandemia afectaron la serie.
Transformación de datos: Si los residuos no son normales, una transformación logarítmica o Box-Cox puede ayudar.
Recordar: El supuesto de normalidad significa que los errores o residuos de un modelo deben seguir una distribución normal (o “campana de Gauss”). Si los errores son normales, podemos hacer predicciones más confiables y usar ciertas pruebas estadísticas que asumen esta propiedad.
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,4)
## Q* = 70.937, df = 18, p-value = 3.135e-08
##
## Model df: 6. Total lags used: 24Paso 4. Pronóstico (modelo manual)
El modelo sigue la tendencia general, pero tiene un sesgo de sobreestimación.
Si la serie tiene patrones estacionales fuertes y no están bien capturados, el modelo puede fallar en prever las fluctuaciones. En ese caso se puede mejorar el modelo, al trabajar desde el inicio con la serie ajustada por estacionalidad en caso de que se detecte un componente estacional fuerte.
Pronóstico en el test de prueba (oct, nov y dic 2024) y gráfico
#Aquí se crea el pronóstico con el modelo ARIMA manual y se guarda en una nueva riable u objeto "manual_forecast"
manual_forecast <- forecast(manual_arima_model, h = length(test_ts)) #Se generan tantos pronósticos como valores tenga el conjunto de prueba (test_ts).
# Crear dataframe para gráfico interactivo del pronóstico manual
manual_forecast_data <- data.frame(Tiempo = time(manual_forecast$mean), ## Extrae las fechas del pronóstico
Pronostico = as.numeric(manual_forecast$mean), ## Valores pronosticados
Observado = as.numeric(test_ts)) ## Valores reales de la serie
# Graficar pronóstico manual junto con los valores observados reales
p_manual <- ggplot(manual_forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico Manual")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Variable1:Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("pncafe")
ggplotly(p_manual)Métricas de evaluación del modelo manual dentro del periodo de prueba (oct,nov y dic2024
# Calcular métricas de evaluación del modelo manual
mae_manual <- mean(abs(manual_forecast$mean - test_ts), na.rm = TRUE)
rmse_manual <- sqrt(mean((manual_forecast$mean - test_ts)^2, na.rm = TRUE))
# Mostrar métricas de evaluación del modelo manual
cat("MAE Manual: ", mae_manual, "\n")## MAE Manual: 547.1456cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 592.5892MAE (Mean Absolute Error)
Indica el error promedio en unidades de la variable, es decir, en número de microempresas.
RMSE (Root Mean Squared Error)
Penaliza más los errores grandes debido a la elevación al cuadrado antes de calcular la raíz.
A continuación se calcula la Tabla de pronóstico modelo manual VS los datos reales u observado en oct,nov y dic2024
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
arima_forecast_manual <- forecast(manual_arima_model, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table_manual <- data.frame(
Tiempo = time(arima_forecast_manual$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(arima_forecast_manual$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table_manual)## Tiempo Observado Pronosticado
## 1 2024.750 1339.155 1108.279
## 2 2024.833 1761.410 1107.805
## 3 2024.917 1798.231 1041.274Ahora pronosticamos fuera del periodo de análisis: Enero 2025
Es decir, le sumamos al periodo de prueba (oct,nov,dic2024) una observación más (enero 2025). Es decir, se estan pronosticando 4 observaciones o meses:
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast_manual <- forecast(manual_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
next_month_forecast_manual <- data.frame(
Tiempo = time(next_forecast_manual$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast_manual$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast_manual)## Tiempo Pronostico
## 1 2024.750 1108.2792
## 2 2024.833 1107.8046
## 3 2024.917 1041.2744
## 4 2025.000 987.3621# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast_manual, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))## [1] "Pronóstico para enero 2025: 2025 = 987.362065627759"Otra forma para calcular un valor futuro (fuera de muestra)-Modelo manual, es decir, en caso de que no se haga la dviisón inicial de conjunto de entrenamiento y prueba
Si no hay conjunto de prueba o test y solo quieres el siguiente punto u observación, el código a usar seria algo asi:
# Pronosticar octubre de 2024
future_forecast_manual <- forecast(manual_arima_model, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024 <- future_forecast_manual$mean[1]
print(paste("Pronóstico para oct 2024:", forecast_oct2024))## [1] "Pronóstico para oct 2024: 1108.27922781123"Modelo ARIMA automático
Usa la función auto.arima() de forecast en R para seleccionar automáticamente los mejores parámetros (p,d,q).
✅ Ventajas:
✔ Optimización automática: Encuentra los valores óptimos de ARIMA sin intervención manual. ✔ Ahorra tiempo: Útil cuando hay muchas series a modelar. ✔ Evita sesgo humano: Reduce el riesgo de elegir un modelo incorrecto por falta de experiencia. ✔ Incluye corrección por estacionalidad si se usa con seasonal = TRUE. ✔ Suele funcionar bien en la mayoría de los casos, ya que usa criterios como AIC/BIC para optimizar.
❌ Desventajas: ❌ Puede no ser el mejor modelo posible, ya que depende del criterio de selección. ❌ Menos interpretabilidad: No siempre es claro por qué eligió ciertos parámetros. ❌ Puede ignorar conocimiento experto sobre la serie o factores externos.
Identificación automática del modelo ARIMA
El modelo automático identificado es (4,1,4). Si se compara el AIC o BIC de este modelo frente el modelo manual (1,1,1), se obtiene un valor más bajo en esta métricas en este modelo automático. Probablemente pudiera ser un buen modelo para pronosticar la variable 1=energia.
library(forecast)
# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)
# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)## Series: train_ts
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.3025 -0.4979
## s.e. 0.0786 0.0773
##
## sigma^2 = 29051: log likelihood = -996.22
## AIC=1998.44 AICc=1998.6 BIC=2007.51
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 14.79646 168.763 133.6028 -0.4502639 12.9936 0.797852 -0.02653252Significancia de coeficientes
library(lmtest)
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.302466 0.078576 -3.8494 0.0001184 ***
## ma2 -0.497857 0.077269 -6.4432 1.17e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Ajuste del modelo ARIMA(4,1,4) 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, 1, 4)) # Especificamos directamente (p=4, d=1, q=4)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(2,1,4)
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4
## -0.5123 -0.6962 0.2398 0.2108 -0.4348 -0.5585
## s.e. 0.1542 0.0922 0.1358 0.0839 0.0714 0.0655
##
## sigma^2 = 26861: log likelihood = -988.64
## AIC=1991.28 AICc=1992.06 BIC=2012.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 14.32523 160.0987 124.9687 -0.3288151 12.23095 0.7462909
## ACF1
## Training set -0.02289033# 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,1,4)
## Q* = 70.937, df = 18, p-value = 3.135e-08
##
## Model df: 6. Total lags used: 24Pronóstico modelo ARIMA automático (2,1,4)
# 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("pncafe")
ggplotly(p4auto) # Convertir el gráfico en interactivoEl modelo automático (2,1,4) 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.
# 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 1339.155 1069.374
## 2 2024.833 1761.410 1031.733
## 3 2024.917 1798.231 1031.733Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueb auna observación más. Es decir, se estan pronosticando 4 observaciones o trimestres.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente trimestre (1 período adicional)
next_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo trimestre
next_month_forecast_auto <- data.frame(
Tiempo = time(next_forecast_auto$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast_auto$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast_auto)## Tiempo Pronostico
## 1 2024.750 1069.374
## 2 2024.833 1031.733
## 3 2024.917 1031.733
## 4 2025.000 1031.733# 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 = 1031.73275916496"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar el mes octubre 2024
future_forecast_auto <- forecast(auto_arima_model_no_seasonal, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024_auto <- future_forecast_auto$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_auto))## [1] "Pronóstico para octubre 2024: 1069.37383824435"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.
Un SARIMA (Seasonal ARIMA) es un modelo que combina:
✅ Un modelo ARIMA tradicional para capturar relaciones entre valores pasados y errores. ✅ Un componente estacional, útil cuando los datos muestran patrones repetitivos en el tiempo (por ejemplo, ventas trimestrales, datos mensuales de temperatura, etc.).
Básicamente, es como un ARIMA mejorado que puede manejar ciclos estacionales.
En este caso, el modelo detectó un patrón repetitivo cada 4 períodos, por lo que usa componentes autoregresivos y de media móvil para modelarlo
El modelo ajustado es un SARIMA(0,1,1)(1,0,0)[12], lo que significa:
(0,1,1): Parte ARIMA no estacional: 0 términos autorregresivos (AR). 1 diferenciación (d), lo que indica que la serie fue diferenciada una vez para hacerla estacionaria. 1 término de media móvil (MA).
(1,0,0)[12]: Parte estacional con periodicidad 12 (mensual si los datos son mensuales): 1 término autorregresivo estacional (SAR). 0 diferenciaciones estacionales. 0 términos de media móvil estacionales (SMA).
El modelo SARIMA(0,1,1)(1,0,0)[12] sugiere que:
- La serie tiene una tendencia no estacionaria, corregida con una diferenciación.
- Existe una influencia significativa del error pasado (MA(1)).
- Hay un componente estacional autorregresivo fuerte cada 12 períodos.
- El ajuste es adecuado según los criterios AIC y BIC, pero se podría comparar con otros modelos para mejorar la predicción.
# Identificación automática del mejor modelo ARIMA
auto_arima_model <- auto.arima(train_ts) # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)## Series: train_ts
## ARIMA(0,1,2)(0,0,2)[12]
##
## Coefficients:
## ma1 ma2 sma1 sma2
## -0.3925 -0.4243 0.3339 0.2808
## s.e. 0.0841 0.0821 0.0872 0.0688
##
## sigma^2 = 23388: log likelihood = -980.05
## AIC=1970.1 AICc=1970.51 BIC=1985.22A continuación, se crea el objeto darima para lueg poder graficar los valores reales y observados:
# Cargar el paquete necesario
library(forecast)
# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12]
darima <- Arima(train_ts,
order = c(2, 1, 4), # (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(2,1,4)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4 sar1
## -0.5118 -0.6617 0.1172 0.1664 -0.3931 -0.5092 0.4531
## s.e. 0.1350 0.1395 0.1286 0.1627 0.0858 0.0779 0.0789
##
## sigma^2 = 22160: log likelihood = -974.79
## AIC=1965.58 AICc=1966.58 BIC=1989.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 11.36139 144.9173 111.9949 -0.267175 10.96431 0.6688137
## ACF1
## Training set -0.05094166En 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(2,1,4)(1,0,0)[12]
## Q* = 32.979, df = 17, p-value = 0.01134
##
## Model df: 7. Total lags used: 24Pronóstico con el modelo SARIMA
# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts)) # Predecir los valores futuros
# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean),
Pronostico = as.numeric(forecast_arima$mean),
Observado = as.numeric(test_ts))
# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
geom_line(aes(y = Pronostico, color = "Pronóstico")) +
geom_line(aes(y = Observado, color = "Observado")) +
ggtitle("Pronóstico vs Observado") +
xlab("Tiempo") + ylab("Unidad Variable 1")
ggplotly(p4) # Convertir el gráfico en 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 1339.155 1186.877
## 2 2024.833 1761.410 1186.487
## 3 2024.917 1798.231 1118.730Ahora pronosticamos fuera del periodo de análisis
Es decir, le sumamos al periodo de prueba una observación más. Es decir, se estan pronosticando 4 observaciones o meses.
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(auto_arima_model, h = length(test_ts) + 1)
# Extraer el pronóstico del próximo mes
next_month_forecast <- data.frame(
Tiempo = time(next_forecast$mean), # Extraer la fecha del pronóstico
Pronostico = as.numeric(next_forecast$mean) # Valor pronosticado
)
# Mostrar el pronóstico completo
print(next_month_forecast)## Tiempo Pronostico
## 1 2024.750 1186.877
## 2 2024.833 1186.487
## 3 2024.917 1118.730
## 4 2025.000 1064.748# 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 = 1064.74767606341"Otra forma para calcular un valor futuro (fuera de muestra)
# Identificación automática del mejor modelo ARIMA
auto_arima_model <- auto.arima(train_ts) # Busca automáticamente los mejores parámetros del modelo ARIMA
print(auto_arima_model)## Series: train_ts
## ARIMA(0,1,2)(0,0,2)[12]
##
## Coefficients:
## ma1 ma2 sma1 sma2
## -0.3925 -0.4243 0.3339 0.2808
## s.e. 0.0841 0.0821 0.0872 0.0688
##
## sigma^2 = 23388: log likelihood = -980.05
## AIC=1970.1 AICc=1970.51 BIC=1985.22# Pronosticar octubre 2024
future_forecast_sarima <- forecast(auto_arima_model, h = 1)
# Extraer el valor específico de octubre 2024
forecast_oct2024_sarima <- future_forecast_sarima$mean[1]
print(paste("Pronóstico para octubre 2024:", forecast_oct2024_sarima))## [1] "Pronóstico para octubre 2024: 1186.87679545544"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 Pablo López & Natalya López