Mayagüez
Antonio Solarte / Juan David Pérez
2025-04-02
library(readxl)
library(tseries)## Warning: package 'tseries' was built under R version 4.4.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)## Warning: package 'forecast' was built under R version 4.4.3library(ggplot2)
library(plotly)## Warning: package 'plotly' was built under R version 4.4.3##
## Adjuntando el paquete: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(timetk)## Warning: package 'timetk' was built under R version 4.4.3library(readxl)
data_col <- read_excel("C:/Users/usuario/Downloads/data_col.xlsx")Introducción
El análisis de series de tiempo es una herramienta fundamental en la toma de decisiones económicas y empresariales, permitiendo anticipar tendencias y mejorar la planificación estratégica. En este estudio, se analizan diversas variables económicas de Cali, como el Índice de Seguimiento a la Economía (ISE) y las exportaciones, para comprender el contexto económico de la región. Sin embargo, las proyecciones se centran exclusivamente en el Índice de Precios al Productor de la Industria Refinadora (IPIR Cali), dado su impacto directo en los costos de producción de empresas como Ingenio Mayagüez, que depende de la estabilidad en los precios de los insumos para la toma de decisiones operativas y financieras. A través de modelos ARIMA, se identificó que el ARIMA(1,0,1) ofrece los pronósticos más precisos, lo que resalta la importancia de aplicar modelos adecuados para obtener estimaciones confiables y útiles para la industria
Variable 1: Índice de seguimiento a la economía nacional (ISE)
variable1_ts <- ts(data_col$ISE, start = c(2012, 1), frequency = 12)
variable1_ts## Jan Feb Mar Apr May Jun Jul
## 2012 81.64043 84.46710 87.79398 84.04757 87.93170 87.83290 88.24695
## 2013 85.11318 86.40213 88.54827 89.80515 92.49118 91.91154 94.40165
## 2014 89.07778 92.67633 95.62235 92.11859 95.55301 95.38251 97.61607
## 2015 91.76237 94.94176 98.53379 95.29803 98.90243 99.26793 102.05825
## 2016 93.39505 99.17002 99.71709 98.04923 100.68981 101.58967 100.45061
## 2017 94.73969 98.80185 102.17733 97.73042 102.05606 104.51569 104.14856
## 2018 96.23811 100.37304 103.83661 101.67362 104.31930 106.35693 107.25107
## 2019 99.73934 103.89290 107.31062 104.10805 108.70366 108.92926 111.83666
## 2020 103.18984 107.69487 100.59924 83.10435 89.76904 95.32219 100.01161
## 2021 99.58701 104.19014 111.84319 104.32340 101.77918 110.07074 113.18207
## 2022 107.74859 111.19639 119.71222 115.06637 118.89908 118.67243 120.19293
## 2023 111.96625 113.81958 121.57124 114.47808 119.31879 121.42211 120.84675
## 2024 113.79232 116.52922 119.14383 120.89703 122.04367 119.72379 125.09974
## Aug Sep Oct Nov Dec
## 2012 87.35799 88.86101 88.77649 93.65018 98.18311
## 2013 93.26500 94.11902 94.27071 99.41349 105.30762
## 2014 97.57139 98.52779 98.45328 102.31746 109.25476
## 2015 101.25747 101.28919 100.69337 104.67015 111.32525
## 2016 104.65810 103.67937 102.02757 107.89716 114.94335
## 2017 105.76854 104.19531 103.26258 108.94782 116.90946
## 2018 109.05707 106.78797 106.89609 112.63571 119.09165
## 2019 111.80893 108.92865 110.32446 115.31801 122.72071
## 2020 100.24538 101.91826 105.08635 110.94001 119.90180
## 2021 111.91996 115.61032 115.79627 122.67239 132.27111
## 2022 121.71072 120.55340 120.06978 122.62911 132.97120
## 2023 122.21000 119.49893 119.21274 126.57906 133.91789
## 2024 124.49061 120.74054 123.58018 127.14149 137.86509Variable 2: Exportaciones totales Cali (X Cali)
variable2_ts <- ts(data_col$X_CALI, start = c(2012, 1), frequency = 12)
variable2_ts## Jan Feb Mar Apr May Jun
## 2012 2776386323 3334123499 3673264590 3098189849 2932225874 2565482794
## 2013 2244740544 2185524243 2464418188 2755741900 2747052270 2647751155
## 2014 2758488920 2843705861 3323962097 2981585788 3271704225 2395820832
## 2015 3174455742 3691354544 4015377640 3472472816 2961046719 3827978560
## 2016 2632054169 2906689482 3434341521 3076631833 2754362479 3307862817
## 2017 3811876173 4648526303 5227303081 4377640917 4746598111 4160808855
## 2018 3887043533 4644828478 3873547690 4928675549 4184578772 4249056636
## 2019 4228509341 3916711568 4648743472 4717923717 4989707006 3674603653
## 2020 3418154940 3685135737 4263706143 5233898931 3818397225 3626080981
## 2021 4278552966 4686144723 5442994707 5496358264 1765523839 3664920178
## 2022 6305393892 3950073657 6726557968 5678729637 5691718620 3887278181
## 2023 4446894806 4072165861 6028195457 4356217896 5109941861 4207398817
## 2024 3167909452 3646367225 3799156534 3951956406 3384485784 4265307769
## Jul Aug Sep Oct Nov Dec
## 2012 3337828773 3394154902 3268970900 3420138034 2914358834 2479710628
## 2013 3009923607 3029354558 3170768901 2995002333 3705333404 3253596203
## 2014 3067058040 3150634166 3607121900 3599117333 3354623853 4149723028
## 2015 3493366175 4544635548 3888813999 3984680418 3914556125 4012245071
## 2016 2349860588 3524679872 5275026533 6480140184 8154599544 6149933370
## 2017 5285983520 4569108625 5338659521 5224985907 5223869913 4683183460
## 2018 3961801153 4466636820 4362039836 5549749086 4967662946 4775473366
## 2019 3861389769 4420553014 4612539042 4429584755 3961324974 4433742707
## 2020 4435913270 4552397354 5186652718 5654196991 4806736880 4549848739
## 2021 5789605073 5888852664 5446599653 3529451909 7540768229 4622189941
## 2022 6346181463 4177626857 6269764101 5445590818 6064310505 7927917277
## 2023 4032178168 4515754300 4520887887 5217843635 4043781413 4054723207
## 2024 4560161323 4043146269 5474405775 4767189703 4664466295 4495265517Variable 3: Índice de producción industrial Cali (IPIR Cali)
variable3_ts <- ts(data_col$IPIR_CALI, start = c(2012, 1), frequency = 12)
variable3_ts## Jan Feb Mar Apr May Jun Jul
## 2012 96.47151 101.89722 107.43327 98.47965 104.51982 106.11454 105.94103
## 2013 93.69983 93.53270 95.12147 106.11457 104.44380 107.49263 111.44949
## 2014 98.13421 101.15981 105.46999 104.37586 104.60483 104.35731 109.90205
## 2015 102.11504 107.40357 111.24739 107.41472 107.90188 110.71140 112.16074
## 2016 104.96732 109.38748 108.97920 107.20438 104.07097 108.11365 103.99047
## 2017 100.22203 102.06038 106.84059 100.92828 103.58373 108.64599 107.97362
## 2018 100.52465 103.07804 107.52431 106.97189 108.88635 108.82576 110.75496
## 2019 104.31989 108.36834 112.78477 110.74097 118.10134 113.66469 114.34179
## 2020 104.42163 110.08244 107.53392 92.19771 98.84423 106.35106 109.87100
## 2021 100.93980 104.76326 107.43094 97.33371 52.56788 101.04241 124.88593
## 2022 105.20084 104.88014 118.25102 113.47250 116.38268 115.92182 121.74192
## 2023 106.59331 108.27840 116.65529 110.15517 117.11739 118.34544 116.78418
## 2024 109.96371 115.47878 114.37051 121.84172 122.02643 119.89085 119.44905
## Aug Sep Oct Nov Dec
## 2012 107.06245 109.51453 116.58799 112.16612 106.41190
## 2013 110.48614 111.73036 116.18060 108.94420 106.24920
## 2014 108.63293 111.80402 116.35420 109.33215 111.67265
## 2015 111.45782 114.56244 118.01970 114.27364 114.38366
## 2016 114.75106 114.25847 115.65555 113.09275 112.05771
## 2017 112.68312 113.10225 116.27296 113.84054 113.11650
## 2018 115.53036 114.93168 119.94612 116.72999 109.58689
## 2019 113.83002 113.34219 117.95999 113.45428 113.20773
## 2020 108.68535 112.61665 116.42163 113.40291 109.39146
## 2021 127.58523 127.32584 124.93936 119.20490 118.04875
## 2022 125.74966 127.96142 129.59093 128.40464 121.89443
## 2023 121.57073 124.15010 124.79659 121.83217 119.30224
## 2024 119.67162 116.32233 118.35843 114.83672 111.36085Tasas de crecimiento 2023-2024
# Cargar librería necesaria
library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Asegurar que la columna Fecha es de tipo Date
data_col$Fecha <- as.Date(data_col$FECHA, format="%Y-%m-%d")
# Filtrar los valores de diciembre 2023 y diciembre 2024
valores_2023 <- data_col %>% filter(format(FECHA, "%Y-%m") == "2023-12")
valores_2024 <- data_col %>% filter(format(FECHA, "%Y-%m") == "2024-12")
# Calcular la tasa de crecimiento para cada variable
tasa_crecimiento <- (valores_2024[, c("IPIR_CALI", "ISE", "X_CALI")] -
valores_2023[, c("IPIR_CALI", "ISE", "X_CALI")]) /
valores_2023[, c("IPIR_CALI", "ISE", "X_CALI")] * 100
# Mostrar la tasa de crecimiento
print(tasa_crecimiento)## IPIR_CALI ISE X_CALI
## 1 -6.656527 2.947475 10.86492Extracción de señales
ISE
Gráfico inicial ISE-Original
library(ggplot2)
library(plotly)data_col$variable1 <- as.numeric(variable1_ts)
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) +
geom_point(color = "black", size = 0.1) +
ggtitle("ISE: Serie original") +
xlab("Tiempo") +
ylab("Puntaje") +
theme_minimal()
ggplotly(grafico_serie)En cuanto a la extracción de señales de la gráfica original del ISE, podemos ver que tiene una tendencia alcista a lo largo del periodo, con una disrupción en la pandemia, la cual frenó en gran medida la economía y por eso podemos ver esa caída tan grande en ese año, sin embargó se recuperó incluso para finales de año, lo que generó que pudiera seguir con su tendencia original. Vemos que tiene una estacionalidad para los último dos meses del año, noviembre y diciembre, y el inicio de año, enero; esto debido a que un gran número de fiestas ocurren en el último trimestre del año, hay mucho más consumo y se dinamiza la economía, sin embargo, los inventarios sobrantes de diciembre quedan estancados para enero y los consumidores ya no tienen la liquidez suficiente para seguir con los volúmenes de compra de meses anteriores. En cuanto al residuo se mantiene estable, durante el periodo de tiempo, lo que indica una ausencia de eventos inesperados para la variable, a excepción de la pandemia.
Extracción señales ISE
library(ggplot2)
library(plotly)stl_decomp_var1 <- stl(variable1_ts, s.window = "periodic")stl_df_var1 <- data.frame(
Time = rep(time(variable1_ts), 4),
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))
)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 ISE",
x = "Tiempo",
y = "Puntaje")ggplotly(p)X Cali
Gráfico inicial X Cali-Original
data_col$variable2 <- as.numeric(variable2_ts)grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month",
length.out = nrow(data_col)),
y = variable2)) +
geom_line(color = "grey", linewidth = 0.4) +
geom_point(color = "black", size = 0.1) +
ggtitle("X Cali: Serie original") +
xlab("Tiempo") +
ylab("Dólares FOB") +
theme_minimal()
ggplotly(grafico_serie)Para el gráfico original de las exportaciones totales de la ciudad de Cali, no se tiene una tendencia bien definida como tal, pero el valor inicial es menor al valor final, con lo que podemos intuir que está subiendo, con unos picos en 2016, 2021, y 2022; y con una caída para 2021, demarcada por el paro nacional de 2021. El residuo tiene un aumento para 2016, cuando hubo un gran crecimiento en las exportaciones de la ciudad, sin embargo, ha estado moviéndose desenfrenadamente sin un patrón visible en la era post-pandemia.
Extracción señales X Cali
stl_decomp_var2 <- stl(variable2_ts, s.window = "periodic")stl_df_var2 <- data.frame(
Time = rep(time(variable2_ts), 4),
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))
)
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 X Cali",
x = "Tiempo",
y = "Dólares FOB")
ggplotly(p)IPIR Cali
Gráfico inicial IPIR Cali-Original
data_col$variable3 <- as.numeric(variable3_ts)grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month",
length.out = nrow(data_col)),
y = variable3)) +
geom_line(color = "grey", linewidth = 0.4) +
geom_point(color = "black", size = 0.1) +
ggtitle("IPIR Cali: Serie original") +
xlab("Tiempo") +
ylab("Puntaje") +
theme_minimal()ggplotly(grafico_serie)El índice de producción industrial regional de Cali se mantuvo estable relativamente entre sus fluctuaciones estacionales para fin de año hasta el paro nacional de 2021, momento donde hay una caída estrepitosa en cuanto a este índice, sin embargo en un par de meses se recuperó y en un nivel superior al que estaba
Extracción señales IPIR Cali
stl_decomp_var3 <- stl(variable3_ts, s.window = "periodic")stl_df_var3 <- data.frame(
Time = rep(time(variable3_ts), 4),
Value = c(stl_decomp_var3$time.series[, "seasonal"],
stl_decomp_var3$time.series[, "trend"],
stl_decomp_var3$time.series[, "remainder"],
variable2_ts),
Component = rep(c("Estacional", "Tendencia", "Residuo", "Serie Original"), each =
length(variable3_ts))
)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 IPIR Cali",
x = "Tiempo",
y = "Puntaje")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:
ISE ajustada por estacionalidad
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]X Cali ajustada por estacionalidad
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]IPIR Cali ajustada por estacionalidad
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]Ahora sí se puede graficar las series originales vs. la ajustada por estacionalidad
# 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## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in geom_line(aes(x = fechas_var1, y = variable1_ts), color = "grey", :
## Ignoring unknown parameters: `name`## Warning in geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", :
## Ignoring unknown parameters: `name`# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.Vemos que cuando ajustamos el gráfico de la serie del ISE con la estacionalidad, tiene un crecimiento constante y creciente hasta la pandemia, y no fue sino hasta mediados de 2021 que se recuperó en términos de ajuste por estacionalidad.
# 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## Warning in geom_line(aes(x = fechas_var2, y = variable2_ts), color = "grey", :
## Ignoring unknown parameters: `name`## Warning in geom_line(aes(x = fechas_var2, y = variable2_sa), color = "black", :
## Ignoring unknown parameters: `name`# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.En la descomposición temporal se mostró una temporalidad para esta variable, sin embargo personalmente no la veo y esto se corroboró con esta gráfica, donde la serie ajustada por estacionalidad es muy idéntica a la original.
# 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## Warning in geom_line(aes(x = fechas_var3, y = variable3_ts), color = "grey", :
## Ignoring unknown parameters: `name`## Warning in geom_line(aes(x = fechas_var3, y = variable3_sa), color = "black", :
## Ignoring unknown parameters: `name`# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.En este gráfico se suavizan las estacionalidades que tiene la producción industrial, y es donde podemos visualizar de mejor medida como la variable estaba en una tendencia horizontal sin cambios drásticos que llamaran la atención, hasta que llegó el paro nacional de 2021, con esa gran caída pero también con esa recuperación un par de meses después, con el aumento de niveles incluso
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
variable1_vec <- as.numeric(variable1_ts)
tendencia_var1 <- as.numeric(stl_decomp_var1$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable1_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var1 <- ggplot() +
geom_line(aes(x = fechas, y = variable1_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var1, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Variable 1: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Unidad de medida Variable 1") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var1)Tendencia Variable 2
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable2_vec <- as.numeric(variable2_ts)
tendencia_var2 <- as.numeric(stl_decomp_var2$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable2_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var2 <- ggplot() +
geom_line(aes(x = fechas, y = variable2_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var2, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Variable 2: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Unidad de medida Variable 2") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var2)Tendencia variable 3
library(ggplot2)
library(plotly)
# Convertir la serie a un vector numérico
variable3_vec <- as.numeric(variable3_ts)
tendencia_var3 <- as.numeric(stl_decomp_var3$time.series[, "trend"])
# Asegurar que 'fechas' tenga la misma longitud
fechas <- seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = length(variable3_ts))
# Gráfico interactivo de la serie original vs tendencia
grafico_tendencia_var3 <- ggplot() +
geom_line(aes(x = fechas, y = variable3_vec, color = "Serie Original"), size = 0.7, linetype = "solid") +
geom_line(aes(x = fechas, y = tendencia_var3, color = "Tendencia"), size = 0.8, linetype = "solid") +
scale_color_manual(values = c("Serie Original" = "grey", "Tendencia" = "black")) +
ggtitle("Variable 3: Serie Original vs Tendencia") +
xlab("Tiempo") +
ylab("Unidad de medida Variable 3") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas del eje X
# Convertir a gráfico interactivo con plotly
ggplotly(grafico_tendencia_var3)Ahora calculamos la tasa de crecimiento de la serie original vs tendencia:
Tasa de crecimiento de la serie de tendencia y original para la variable 1
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var1 <- (variable1_ts[(13:length(variable1_ts))] / variable1_ts[1:(length(variable1_ts) - 12)] - 1) * 100
tasa_tendencia_var1 <- (tendencia_var1[(13:length(tendencia_var1))] / tendencia_var1[1:(length(tendencia_var1) - 12)] - 1) * 100
# Crear vector de fechas corregido, es decir que inicie desde enero 2013
fechas_corregidas_var1 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var1))
# Verificar longitudes
print(length(fechas_corregidas_var1))## [1] 144print(length(tasa_crecimiento_var1))## [1] 144print(length(tasa_tendencia_var1))## [1] 144Gráfico variable original y tendencia variable 1: tasa de crecimiento anual
library(ggplot2)
library(plotly)
# Gráfico de la tasa de crecimiento anual variable 1
grafico_crecimiento_var1 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var1, y = tasa_crecimiento_var1), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable1: Tasa de crecimiento anual % de la serie Original y la tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)Vemos que la tasa de crecimiento anual para la variable ISE venía decreciendo poco a poco y acercándose cada vez más al estancamiento, podemos observar que la pandemia generó que esa tasa se convirtiera en negativa pero hubo una rápida recuperación, en donde se vieron tasas de crecimiento anual por encima del 10% para mediados de 2021, algo no visto antes, sin embargo para mediados de 2023 volvió a sus estado previo donde las tasas son casi estáticas.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))
# Verificar longitudes
print(length(fechas_corregidas_var2))## [1] 144print(length(tasa_crecimiento_var2))## [1] 144print(length(tasa_tendencia_var2))## [1] 144# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var2 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_crecimiento_var2), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var2, y = tasa_tendencia_var2), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable2: Tasa de crecimiento anual % de la serie Original y la Tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)Para las exportaciones totales de Cali, hay una preocupante reincidencia en tasas de crecimiento anual negativas, momentos donde incluso se han alcanzado los -20%, algo malo para el sector exportador de la ciudad. También se han visto algunos periodos prósperos para el sector donde las tasas de crecimiento han alcanzado cerca del 50%.
Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 3
#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100
# Crear vector de fechas corregido
fechas_corregidas_var3 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var3))
# Verificar longitudes
print(length(fechas_corregidas_var3))## [1] 144print(length(tasa_crecimiento_var3))## [1] 144print(length(tasa_tendencia_var3))## [1] 144# Gráfico de la tasa de crecimiento anual variable 2
grafico_crecimiento_var3 <- ggplot() +
geom_line(aes(x = fechas_corregidas_var3, y = tasa_crecimiento_var3), color = "grey", size = 0.7) +
geom_line(aes(x = fechas_corregidas_var3, y = tasa_tendencia_var3), color = "black", size = 0.8, linetype = "dashed") +
ggtitle("Variable3: Tasa de crecimiento anual % de la serie Original y la tendencia") +
xlab("Tiempo") +
ylab("% de Crecimiento Anual") +
theme_minimal()
# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)Y finalmente, la tasa de crecimiento de la producción industrial regional de la ciudad se ha mantenido estancado durante casi todo el periodo a evaluar, con ligeros tasas negativas en pandemia, tasas positivas la era posterior al paro nacional de 2021, y la actualidad, que cuenta con una tasa de crecimiento negativa para la producción industrial.
Modelo ARIMA manual
Identificar estacionariedad
Se escogio la variable 3 (IPIR CALI)
library(tseries)# 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(variable3_ts) - 3 # Se deja fuera los últimos 3 valores para usarlos como set de prueba.
train_ts <- window(variable3_ts, end = c(2024, 9)) # Entrenamiento hasta septiembre 2024
test_ts <- window(variable3_ts, start = c(2024, 10)) # Prueba inicia desde oct2024library(tseries)
# Prueba de estacionariedad con Augmented Dickey-Fuller (ADF)
adf_test <- adf.test(train_ts) #Se aplica el test ADF a la variable 1 (conjunto de entrenamiento)## Warning in adf.test(train_ts): p-value smaller than printed p-valueprint(adf_test) # se muestra el resultado del test##
## Augmented Dickey-Fuller Test
##
## data: train_ts
## Dickey-Fuller = -5.818, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryp-value = 0.01: Como el p-valor es menor a 0.05, rechazamos la hipótesis nula de que la serie tiene una raíz unitaria (es decir, que no es estacionaria).
La serie train_ts es estacionaria, por lo que no es necesario diferenciarla antes de ajustar un modelo ARIMA.
#Se crea un nuevo objeto o variable que se llama train_diff, en donde se diferencia la variable 3 , una sola vez:
train_diff <- diff(train_ts, differences = 1) # Graficar la serie original en un gráfico separado
p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable3 = as.numeric(train_ts)), aes(x = Tiempo, y = variable3)) +
geom_line(color = "blue") +
ggtitle("Variable 3:Serie Original") +
xlab("Tiempo") + ylab("Gwh")
ggplotly(p2) # Convertir en gráfico interactivo## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.Tendencia y Estacionalidad: La serie parece mostrar cierta tendencia creciente con fluctuaciones. También hay posibles patrones cíclicos que podrían indicar estacionalidad.
Choque en 2020: Hay una caída abrupta alrededor de 2020, que podría deberse a un evento externo (por ejemplo, la pandemia de COVID-19 afectando el consumo de energía).
# Graficar la serie diferenciada en un gráfico separado
p3 <- ggplot(data.frame(Tiempo = time(train_ts)[-1], variable3_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable3_Diff)) +
geom_line(color = "red") +
ggtitle("variable 3:Serie Estacionaria (Una diferenciación)") +
xlab("Tiempo") + ylab("variable 3 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivo# Si la serie no es estacionaria (p-valor > 0.05), aplicar diferenciación
train_diff_log <- diff(log(train_ts), differences = 1) # 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("Variable1:Serie Estacionaria (Una diferenciación en logaritmo)") +
xlab("Tiempo") + ylab("Variable 3 diferenciada")
ggplotly(p3) # Convertir en gráfico interactivo# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
adf_test_diff <- adf.test(train_diff)## Warning in adf.test(train_diff): p-value smaller than printed p-value print(adf_test_diff)##
## Augmented Dickey-Fuller Test
##
## data: train_diff
## Dickey-Fuller = -7.6193, 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)## Warning in adf.test(train_diff_log): p-value smaller than printed p-value print(adf_test_diff_log)##
## Augmented Dickey-Fuller Test
##
## data: train_diff_log
## Dickey-Fuller = -7.517, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryIdentificación manual de p y q
library(forecast)
# Graficar ACF y PACF
acf_plot <- ggAcf(train_diff_log, lag.max = 6) + ggtitle("Autocorrelation Function (ACF)-Determinar q")
pacf_plot <- ggPacf(train_diff_log, lag.max = 6) + ggtitle("Partial Autocorrelation Function (PACF)-Determinar p")
ggplotly(acf_plot)La autocorrelación cae rápidamente y los valores están dentro del intervalo de confianza.Esto indica que no hay una estructura clara de cola larga, lo que sugiere que un modelo MA de orden bajo podría ser apropiado.
El primer rezago es significativo, por lo que un MA(1) podría ser una buena opción
ggplotly(pacf_plot)Interpretación correlogramas
Se puede observar que los valores que podrian tomar p y q serian:
p=1* p=2* q=1* q=2 (P optimo=1,2) (q óptimo=1)
El modelo óptimo para esta variable seria (1,0,1)
Estimación manual del modelo
# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(1,0,1)) #Se va a estimar un modelo Arima de orden (1,0,1)
summary(manual_arima_model)## Series: train_ts
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.5113 0.1855 110.6325
## s.e. 0.1115 0.1282 1.3372
##
## sigma^2 = 48.14: log likelihood = -512.22
## AIC=1032.44 AICc=1032.71 BIC=1044.56
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.06329217 6.870193 4.341643 -0.4971302 4.283514 0.8828345
## ACF1
## Training set -0.007609028MAPE < 10% → Muy buen modelo
Significancia de coefientes
library(lmtest)## Warning: package 'lmtest' was built under R version 4.4.3## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# 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.51134 0.11152 4.5852 4.536e-06 ***
## ma1 0.18549 0.12825 1.4464 0.1481
## intercept 110.63248 1.33716 82.7370 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El modelo ARIMA(1,0,1) muestra que la serie tiene una fuerte dependencia con su valor inmediato anterior (AR(1) = 0.5113, significativo), mientras que el efecto de los choques pasados es débil (MA(1) = 0.1855, no significativo). La serie se estabiliza alrededor de una media de 110.63 y presenta un error medio absoluto de 4.34 unidades, con un MAPE del 4.28%, indicando un buen ajuste. Además, los residuos no muestran autocorrelación, lo que sugiere que el modelo captura bien la estructura de la serie.
Validación de residuos del modelo estimado manual
checkresiduals(manual_arima_model)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 43.259, df = 22, p-value = 0.004385
##
## Model df: 2. Total lags used: 241. Serie de residuos (gráfico superior)
Los residuos parecen oscilar alrededor de cero sin una tendencia clara, lo que sugiere que el modelo ha capturado bien la estructura de la serie. Sin embargo, hay un pico fuerte alrededor de 2020, lo que podría indicar un evento atípico o un cambio estructural en los datos
Función de Autocorrelación (gráfico inferior izquierdo, ACF de residuos)
La mayoría de los valores están dentro de las bandas de confianza (líneas azules), lo que indica que no hay una autocorrelación significativa en los residuos. Esto sugiere que el modelo ha logrado capturar la dependencia temporal correctamente. Sin embargo, hay algunos picos que podrían indicar pequeños patrones no modelados
Histograma de residuos con ajuste normal (gráfico inferior derecho)
La distribución de los residuos es aproximadamente normal, pero con colas ligeramente más pesadas de lo esperado (algunos valores extremos), lo que sugiere la presencia de valores atípicos. Esto podría afectar la precisión de algunas predicciones
Pronóstico (modelo manual)
#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("Variable3:Pronóstico Manual vs Observado") +
xlab("Tiempo") + ylab("Variable3")
ggplotly(p_manual)## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.La gráfica compara el pronóstico manual del modelo ARIMA(1,0,1) con los valores observados de la serie temporal para la Variable 3, mostrando una tendencia descendente en ambos casos. Sin embargo, el modelo tiende a subestimar los valores reales, ya que el pronóstico es sistemáticamente más bajo que lo observado. A pesar de esto, en el último punto ambas líneas convergen, lo que sugiere que el modelo podría ajustarse mejor en el corto plazo. Para mejorar la precisión, sería recomendable ajustar los parámetros o explorar otras configuraciones del modelo
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: 2.494705cat("RMSE Manual: ", rmse_manual, "\n")## RMSE Manual: 3.150609Se presenta un MAE de 2.49, lo que indica que, en promedio, las predicciones se desvían en aproximadamente 2.49 unidades del valor real. El RMSE de 3.15 sugiere que los errores individuales pueden ser mayores en algunos casos, ya que penaliza más los errores grandes. Aunque estos valores indican un desempeño relativamente aceptable, la diferencia entre MAE y RMSE sugiere que existen algunas desviaciones más grandes en ciertas predicciones
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 118.3584 113.6069
## 2 2024.833 114.8367 112.1534
## 3 2024.917 111.3609 111.4102Se observa que el valor real disminuye progresivamente de 118.35 a 111.36, mientras que el modelo también predice una tendencia descendente, pero con valores más bajos que los observados en los primeros periodos. Esto indica que el modelo capta la dirección general de la serie, aunque subestima los valores en algunos puntos. La diferencia entre el valor observado y el pronosticado refleja el error del modelo, lo que se debe tener en cuenta para evaluar su precisión y posibles ajustes
Ahora 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 113.6069
## 2 2024.833 112.1534
## 3 2024.917 111.4102
## 4 2025.000 111.0302# 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 = 111.030159683933"Modelo ARIMA automático
library(forecast)
# Ajustar un modelo ARIMA automático sin estacionalidad, por eso se pone seasonal=FALSE
auto_arima_model_no_seasonal <- auto.arima(train_ts, seasonal = FALSE)
# Mostrar el modelo seleccionado
summary(auto_arima_model_no_seasonal)## Series: train_ts
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## 0.6311 -0.1583 -0.9461
## s.e. 0.0836 0.0822 0.0299
##
## sigma^2 = 47.48: log likelihood = -508.26
## AIC=1024.52 AICc=1024.79 BIC=1036.62
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.9352405 6.799612 4.343458 0.3336447 4.256151 0.8832036
## ACF1
## Training set -0.02255257Tiene un RMSE de 6.80 y un MAE de 4.34, que son más altos que los obtenidos con el modelo manual (RMSE = 3.15, MAE = 2.49). Esto indica que el modelo manual tuvo un menor error en las predicciones, lo que sugiere que pudo haber captado mejor la estructura de la serie de tiempo en este caso específico
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.631109 0.083610 7.5482 4.412e-14 ***
## ar2 -0.158307 0.082164 -1.9267 0.05401 .
## ma1 -0.946054 0.029892 -31.6493 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Ajuste del modelo ARIMA(4,1,2) automático sin parte estacional y crearlo como variable darima_auto para luego poder graficarlo y crear la tabla
darima_auto <- Arima(train_ts,
order = c(4, 0, 2)) # Especificamos directamente (p=4, d=1, q=2)
# Mostrar resumen del modelo ajustado
summary(darima_auto)## Series: train_ts
## ARIMA(4,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 mean
## 1.0776 -1.2536 0.7282 -0.1526 -0.3819 0.8489 110.6425
## s.e. 0.2619 0.2393 0.1180 0.0951 0.2584 0.1206 1.3292
##
## sigma^2 = 48.15: log likelihood = -510.27
## AIC=1036.54 AICc=1037.54 BIC=1060.78
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0581532 6.778513 4.310128 -0.4899056 4.244362 0.8764262
## ACF1
## Training set -0.0132485El modelo ARIMA(4,1,2) muestra que los coeficientes ar1 y ma1 son los más significativos, con valores de 1.3739 y -1.7813 respectivamente, lo que indica una fuerte dependencia de la serie con sus valores pasados
# 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(4,0,2) with non-zero mean
## Q* = 33.431, df = 18, p-value = 0.01479
##
## Model df: 6. Total lags used: 24La serie de residuos a lo largo del tiempo muestra una distribución relativamente estable, salvo por una perturbación importante alrededor de 2020, lo que podría indicar un evento atípico o una ruptura estructural en la serie. En el gráfico de autocorrelación (ACF), la mayoría de los rezagos se encuentran dentro del intervalo de confianza, aunque algunos picos aislados superan los límites, lo cual sugiere que persiste cierta autocorrelación en los residuos y que el modelo podría no haber captado completamente toda la estructura temporal. Finalmente, el histograma de los residuos con la curva de densidad superpuesta indica una distribución aproximadamente normal, aunque con colas pesadas
Pronóstico modelo ARIMA automático (4,1,2)
# 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("variable3")
ggplotly(p4auto) # Convertir el gráfico en interactivo## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.Ambos presentan una tendencia descendente, lo que indica que el modelo logró identificar correctamente la dirección general del comportamiento. Sin embargo, se observa una discrepancia en la magnitud de la caída: los valores observados descienden de manera más pronunciada, mientras que el modelo suaviza esa caída. En particular, al final del periodo, el pronóstico resulta más alto que el valor real, lo que indica que el modelo subestimó la caída final. Esto sugiere que, aunque el ARIMA logró ajustarse en términos de tendencia, no fue lo suficientemente sensible para captar la intensidad del descenso reciente, posiblemente debido a una estructura del modelo que suaviza excesivamente los cambios bruscos en la serie
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
forecast_arima_auto <- forecast(darima_auto, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_arima_auto <- data.frame(
Tiempo = time(forecast_arima_auto$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(forecast_arima_auto$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_arima_auto)## Tiempo Observado Pronosticado
## 1 2024.750 118.3584 113.3648
## 2 2024.833 114.8367 112.5025
## 3 2024.917 111.3609 111.9925Ahora 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(darima_auto, 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 113.3648
## 2 2024.833 112.5025
## 3 2024.917 111.9925
## 4 2025.000 110.8812# 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 = 110.881168607505"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar el mes octubre 2024
future_forecast_auto <- forecast(darima_auto, 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: 113.364757255794"Modelo SARIMA automático
# 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,3)(0,0,1)[12]
##
## Coefficients:
## ma1 ma2 ma3 sma1
## -0.4136 -0.3406 -0.1595 0.3035
## s.e. 0.0824 0.0861 0.0774 0.0791
##
## sigma^2 = 43.59: log likelihood = -501.76
## AIC=1013.51 AICc=1013.93 BIC=1028.63# Cargar el paquete necesario
library(forecast)
# Ajustar el modelo SARIMA(0,1,1)(1,0,0)[12]
darima <- Arima(train_ts,
order = c(1, 0, 1), # (p,d,q) -> (1,0,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(1,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sar1 mean
## 0.4684 0.1104 0.3362 110.7584
## s.e. 0.1325 0.1500 0.0805 1.5806
##
## sigma^2 = 43.27: log likelihood = -504.18
## AIC=1018.36 AICc=1018.76 BIC=1033.51
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.1246673 6.491285 3.823062 -0.4223033 3.809749 0.7773857
## ACF1
## Training set -0.006425196# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima) # Verificar si los residuos son aleatorios y no presentan patrones
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1)(1,0,0)[12] with non-zero mean
## Q* = 22.587, df = 21, p-value = 0.3664
##
## Model df: 3. Total lags used: 24En la parte superior, la serie de residuos muestra una distribución centrada en cero a lo largo del tiempo, con excepción de un pico abrupto alrededor de 2020, probablemente asociado a un evento atípico o choque externo. El gráfico de autocorrelación (ACF) de los residuos revela que la mayoría de los rezagos están dentro del intervalo de confianza, aunque algunos puntos aislados superan ligeramente los límites, lo cual sugiere que persiste una leve autocorrelación no explicada por el modelo. Por su parte, el histograma muestra que los residuos siguen una distribución aproximadamente normal, aunque con ligera asimetría y presencia de valores extremos en ambos extremos.
Pronó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## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.Ambos presentan una tendencia descendente, pero el pronóstico muestra una caída menos pronunciada que la observada en los datos reales. Esto sugiere que el modelo logra capturar la dirección general de la serie, pero subestima la magnitud del descenso reciente, lo que puede indicar que no está reaccionando completamente a cambios más abruptos o recientes en la dinámica de la serie. A pesar de esta diferencia, la cercanía entre ambas líneas indica un ajuste razonable, aunque podría mejorarse incluyendo más sensibilidad a variaciones de corto plazo o reforzando el componente estacional si este sigue teniendo relevancia. En general, el modelo ofrece un pronóstico conservador frente a una realidad que presenta una caída más acelerada
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Generar pronóstico con el modelo ARIMA identificado
forecast_arima <- forecast(darima, h = length(test_ts))
# Crear un dataframe con los valores observados y pronosticados
forecast_table <- data.frame(
Tiempo = time(forecast_arima$mean), # Extraer las fechas del pronóstico
Observado = as.numeric(test_ts), # Valores reales
Pronosticado = as.numeric(forecast_arima$mean) # Valores pronosticados
)
# Mostrar la tabla
print(forecast_table)## Tiempo Observado Pronosticado
## 1 2024.750 118.3584 115.7977
## 2 2024.833 114.8367 114.6310
## 3 2024.917 111.3609 113.7008Ahora pronosticamos fuera del periodo de análisis
# Cargar librerías necesarias
library(forecast)
# Hacer un pronóstico para el siguiente mes (1 período adicional)
next_forecast <- forecast(darima, 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 115.7977
## 2 2024.833 114.6310
## 3 2024.917 113.7008
## 4 2025.000 110.5241# 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 = 110.524097022942"Otra forma para calcular un valor futuro (fuera de muestra)
# Pronosticar octubre 2024
future_forecast_sarima <- forecast(darima, 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: 115.797684306107"Relación entre el IPIR Cali y la empresa Ingenio Mayagüez
El Índice de Precios al Productor Industrial Regional (IPIR) de Cali es un indicador clave para Ingenio Mayagüez, ya que refleja la evolución de los costos de producción en la región donde opera la empresa. A través del análisis de modelos de series temporales, se encontró que el modelo ARIMA (1,0,1) manual fue el más preciso para predecir el comportamiento del IPIR Cali. Este modelo le permite a la empresa prever tendencias en los costos de insumos esenciales como fertilizantes, combustibles y mano de obra, lo que facilita la toma de decisiones estratégicas en la planeación financiera y operativa. Si el modelo proyecta un aumento en el IPIR Cali, Ingenio Mayagüez podría anticipar incrementos en sus costos de producción y ajustar su estrategia de compra de insumos o precios de venta para proteger su rentabilidad. Por el contrario, si se pronostica una disminución en el índice, la empresa podría aprovechar mejores condiciones de mercado para optimizar su estructura de costos. La precisión del modelo manual ARIMA (1,0,1) refuerza la importancia de un análisis detallado y personalizado en la predicción de indicadores económicos clave, brindándole a Mayagüez una ventaja competitiva en la planificación y gestión de su negocio.
Conclusión
El uso del modelo ARIMA (1,0,1) para analizar y predecir el IPIR Cali demuestra la utilidad de las herramientas de modelado estadístico en la toma de decisiones empresariales. Para Ingenio Mayagüez, este enfoque representa una oportunidad de anticiparse a cambios en los costos de producción, optimizar su cadena de suministro y ajustar su estrategia de precios con base en datos precisos. La elección del modelo manual sobre el automático resalta la importancia de adaptar los métodos analíticos a las características específicas de los datos, garantizando así una mayor precisión en las predicciones. En un entorno de alta volatilidad económica, contar con modelos confiables como este permite a las empresas agroindustriales fortalecer su resiliencia y mantener su competitividad en el mercado.