Series de tiempo: Extracción de señales y Modelo ARIMA

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)
data_col <- read_excel("Base Caso2 (1).xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))
View(data_col)

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

Variable 1

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

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2012 1595.273 1544.653 1718.766 1574.681 1687.182 1630.558 1673.152 1700.788
## 2013 1666.811 1527.898 1648.842 1722.114 1745.650 1637.149 1727.935 1731.272
## 2014 1706.919 1677.512 1765.809 1676.925 1783.913 1700.258 1775.461 1773.788
## 2015 1722.385 1669.782 1805.358 1672.252 1824.030 1719.355 1799.682 1779.533
## 2016 1766.639 1752.126 1747.462 1685.912 1745.006 1712.774 1686.072 1823.924
## 2017 1720.523 1653.542 1747.670 1667.423 1748.151 1706.837 1736.774 1792.472
## 2018 1774.291 1670.436 1800.131 1769.603 1828.144 1789.760 1850.545 1922.534
## 2019 1835.277 1739.269 1884.378 1816.874 1916.032 1836.890 1904.306 1953.351
## 2020 1910.636 1886.473 1793.435 1370.977 1562.429 1616.347 1807.862 1827.641
## 2021 1772.665 1724.037 1940.417 1847.864 1730.073 1871.616 2001.046 2007.131
## 2022 2015.024 1910.252 2105.817 2015.557 2109.814 2049.652 2137.146 2175.928
## 2023 2111.918 1966.525 2195.879 2039.239 2162.926 2076.371 2109.956 2219.682
## 2024 2111.494 2035.827 2120.471 2080.841 2116.916 2008.870 2114.820 2110.274
## 2025 2110.975 2018.777 2174.954 2035.961 2192.316 2118.460 2249.288 2234.768
##           Sep      Oct      Nov      Dec
## 2012 1701.881 1702.813 1667.827 1602.992
## 2013 1697.606 1749.571 1715.983 1666.543
## 2014 1768.670 1809.629 1713.567 1711.971
## 2015 1827.011 1849.687 1736.973 1780.939
## 2016 1739.739 1757.674 1716.030 1722.647
## 2017 1748.235 1776.851 1741.672 1719.700
## 2018 1836.883 1885.882 1838.045 1833.789
## 2019 1884.420 1935.010 1860.997 1926.223
## 2020 1800.417 1853.505 1790.791 1816.589
## 2021 2022.873 2094.568 2032.648 2013.006
## 2022 2117.970 2178.745 2088.379 2073.252
## 2023 2172.905 2191.480 2118.676 2096.022
## 2024 2080.620 2155.514 2061.612 2085.775
## 2025 2207.604 2266.162 2213.671 2218.968

Variable 2

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

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012  81.68037  84.51632  87.83168  84.06845  87.93712  87.82775  88.25721
## 2013  85.15456  86.44792  88.57971  89.82525  92.49853  91.90832  94.41630
## 2014  89.11942  92.72376  95.65359  92.13539  95.55825  95.38099  97.62941
## 2015  91.80939  94.99114  98.56449  95.31551  98.91123  99.26498 102.07327
## 2016  93.44710  99.21948  99.74312  98.06612 100.70060 101.58869 100.46836
## 2017  94.79901  98.86261 102.21842  97.75056 102.06451 104.51566 104.16710
## 2018  96.28252 100.42103 103.86803 101.68526 104.31974 106.34993 107.25595
## 2019  99.80323 103.95610 107.34393 104.13618 108.71167 108.95077 111.85249
## 2020 103.18987 107.70887 100.60622  83.11103  89.76884  95.29142  99.95936
## 2021  99.70037 104.28788 111.86952 104.27038 101.68314 109.96299 113.15155
## 2022 107.67326 111.27952 119.70560 115.14803 118.72427 118.50594 119.94641
## 2023 112.30882 114.01867 121.40610 114.27288 119.03797 121.46259 120.67058
## 2024 114.16295 116.73418 118.87851 120.61068 121.59824 119.06014 124.33195
## 2025 116.39258 118.21867 124.37117 120.82200 124.56234 123.37864 130.23594
##            Aug       Sep       Oct       Nov       Dec
## 2012  87.37441  88.85174  88.73231  93.59978  98.11226
## 2013  93.28462  94.10808  94.22677  99.36308 105.23580
## 2014  97.58678  98.51913  98.41323 102.27108 109.18026
## 2015 101.27515 101.27730 100.65431 104.62138 111.24185
## 2016 104.68155 103.66200 101.98883 107.84345 114.85774
## 2017 105.78888 104.17066 103.21316 108.88465 116.81808
## 2018 109.07291 106.77434 106.87439 112.60466 119.00842
## 2019 111.83213 108.91560 110.26942 115.24839 122.60133
## 2020 100.21522 101.87361 105.09800 111.01645 119.94404
## 2021 111.92693 115.57074 115.80542 122.71004 132.30682
## 2022 121.65503 120.45365 120.52080 122.80948 133.00024
## 2023 122.43384 119.86254 119.71535 127.43698 134.24046
## 2024 124.03954 120.49445 123.64984 127.65571 138.80816
## 2025 127.13485 125.09699 126.95107 130.71397 141.24433

Variable 3

# Convertir/declarar las exportaciones de combustibles en serie de tiempo mensual
variable3_ts <- ts(data_col$X_COMB, start = c(2012, 1), frequency = 12)
variable3_ts

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2012 3313842.2 3319381.0 3870610.8 3522192.3 3529123.0 2884610.0 3198849.0
## 2013 3242512.0 3008895.1 3015949.2 3095551.7 3494107.5 3263329.5 3092261.1
## 2014 3514767.7 2777860.3 2892045.0 2795103.8 3776531.4 3219858.8 3295817.8
## 2015 1564542.2 1787638.6 2008364.4 1769222.8 1948457.0 1779805.3 1539580.2
## 2016  892320.2  998451.8 1006425.5 1050889.4 1378526.2 1481107.8 1303436.7
## 2017 1606633.2 1417686.8 1660507.6 1463076.6 1880608.2 1424209.0 1654446.6
## 2018 1958046.0 1620569.0 1895704.7 2201645.8 2076817.5 1874732.7 2177571.0
## 2019 1615647.9 1724779.6 1915446.7 2305057.2 2200760.2 1817293.5 1744232.7
## 2020 2043776.8 1483145.6  950857.2  701773.5  928404.8  921343.9  936559.7
## 2021 1182946.7 1252154.3 1399605.3 1203308.6 1549776.4 1558797.1 1415764.2
## 2022 2044351.9 1992883.2 2452788.8 3271666.7 2398280.9 3221138.1 3649554.3
## 2023 1986087.5 2314800.3 2185496.0 1834740.4 2107092.1 1966557.4 2109701.2
## 2024 1758565.4 1767110.1 1880150.3 2053557.1 2065070.9 1786374.4 2227329.9
## 2025 1479980.5 1374729.2 1704501.5 1344704.1 1617917.1 1468760.4 1587043.3
##            Aug       Sep       Oct       Nov       Dec
## 2012 2728452.1 3207680.9 3594652.8 2998226.8 3301862.5
## 2013 3523905.9 3296907.8 3208409.6 3382788.8 3658324.0
## 2014 3222883.3 3400174.9 2706862.8 2237362.9 2160977.6
## 2015 1521449.5 1415501.8 1476355.5 1200279.2 1176466.2
## 2016 1542786.3 1373918.1 1395529.2 1333395.9 1687106.1
## 2017 1639329.4 1894853.6 1800861.2 1664315.6 2596000.5
## 2018 2136911.3 2145631.9 2270315.4 1973431.1 1966759.1
## 2019 1704014.1 1616732.8 1712817.4 1526200.8 1748621.8
## 2020 1103433.3  931901.7  969246.6  872364.4 1092089.4
## 2021 1513093.1 1854123.2 1890737.8 2062251.5 2297288.1
## 2022 2225049.8 2590601.0 2183740.1 2544083.4 2456611.3
## 2023 1921370.0 2176384.6 2221638.3 2087084.6 2337004.3
## 2024 1656952.7 1791341.1 1767167.6 1849219.8 1857556.4
## 2025 1214240.6 1593708.2 1389501.2 1316741.9 1430491.5

Extracció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 creación de microempresas en Cali, 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 una serie financiera, 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 economía.

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 en niveles -Original

library(ggplot2)
library(plotly)

# Convertir la serie temporal a un vector numérico para lograr graficar con ggplot2
data_col$variable1 <- as.numeric(variable1_ts)

# Crear el gráfico
grafico_serie <- ggplot(data_col, aes(x = seq.Date(from = as.Date("2012-01-01"), by = "month", length.out = nrow(data_col)), 
                                      y = variable1)) +
  geom_line(color = "grey", linewidth = 0.4) +  # Cambiado 'size' por 'linewidth'
  geom_point(color = "black", size = 0.1) +
  ggtitle("Variable 1: Serie original") +
  xlab("Tiempo") +
  ylab("Unidad Variable 1") +
  theme_minimal()

ggplotly(grafico_serie)

Extracción señales variable 1

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

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

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

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

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

Extracción señales variable 2

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

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

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

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

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

Extracción señales variable 3

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

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

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

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

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

Después de la descomposición temporal de cada variable, se extrae la variable ajustada por estacionalidad para graficarla junto con la serie original:

Se crea la variable1 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$time.series[, "seasonal"]

Se crea la variable2 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable2_sa <- variable2_ts - stl_decomp_var2$time.series[, "seasonal"]

Se crea la variable3 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable3_sa <- variable3_ts - stl_decomp_var3$time.series[, "seasonal"]

Ahora si se puede graficar las series originales versus la ajustada por estacionalidad

Gráfico serie original VS ajustada Variable 1

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

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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)

Gráfico serie original VS ajustada Variable 2

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

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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var2)

Gráfico serie original VS ajustada Variable 3

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

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

# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var3)

Ahora graficamos serie original vs tendencia

• 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] 156

print(length(tasa_crecimiento_var1))

## [1] 156

print(length(tasa_tendencia_var1))

## [1] 156

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

library(ggplot2)
library(plotly)

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

Ahora calculamos la tasa de crecimiento de la serie original vs tendencia: variable 2

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var2 <- (variable2_ts[(13:length(variable2_ts))] / variable2_ts[1:(length(variable2_ts) - 12)] - 1) * 100
tasa_tendencia_var2 <- (tendencia_var2[(13:length(tendencia_var2))] / tendencia_var2[1:(length(tendencia_var2) - 12)] - 1) * 100

# Crear vector de fechas corregido
fechas_corregidas_var2 <- seq(from = as.Date("2013-01-01"), by = "month", length.out = length(tasa_crecimiento_var2))

# Verificar longitudes
print(length(fechas_corregidas_var2))

## [1] 156

print(length(tasa_crecimiento_var2))

## [1] 156

print(length(tasa_tendencia_var2))

## [1] 156

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

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] 156

print(length(tasa_crecimiento_var3))

## [1] 156

print(length(tasa_tendencia_var3))

## [1] 156

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

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

1. Autoregressive o autorregresivo (AR) → Usa valores pasados para predecir el futuro
2. 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.

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

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

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

Paso 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 variable 1 (conjunto de entrenamiento)
print(adf_test) # se muestra el resultado del test

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_ts
## Dickey-Fuller = -3.061, Lag order = 5, p-value = 0.1342
## alternative hypothesis: stationary

El 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:

 # Graficar la serie original en un gráfico separado
  p2 <- ggplot(data.frame(Tiempo = time(train_ts), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
    geom_line(color = "blue") +
    ggtitle("Variable 1:Serie Original") +
    xlab("Tiempo") + ylab("Gwh")
  
  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], variable1_Diff = as.numeric(train_diff)), aes(x = Tiempo, y = variable1_Diff)) +
    geom_line(color = "red") +
    ggtitle("variable 1:Serie Estacionaria (Una diferenciación)") +
    xlab("Tiempo") + ylab("variable 1 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

Ejemplo 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), variable1 = as.numeric(train_ts)), aes(x = Tiempo, y = variable1)) +
    geom_line(color = "blue") +
    ggtitle("Variable1:Serie Original") +
    xlab("Tiempo") + ylab("Variable1")
  
  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("Variable1:Serie Estacionaria (Una diferenciación en logaritmo)") +
    xlab("Tiempo") + ylab("Variable 1 diferenciada")
  
  ggplotly(p3)  # Convertir en gráfico interactivo

Ahora 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.3555, 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.4411, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

En 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_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)

ggplotly(pacf_plot)

Interpretación correlogramas

Se puede observar que los valores que podrian tomar p y q serian:

p=1* p=4 p=6 q=1* q2 q=5 q=6 (P optimo=1) (q óptimo=1)

El modelo óptimo para esta variable seria (1,1,1)

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 (1,1,1)

# Cálculo manual de modelo ARIMA
manual_arima_model <- Arima(train_ts, order = c(1,1,1)) #Se va a estimar un modelo Arima de orden (1,1,1)
summary(manual_arima_model)

## Series: train_ts 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.4590  -0.0642
## s.e.   0.1695   0.2002
## 
## sigma^2 = 5805:  log likelihood = -873.48
## AIC=1752.95   AICc=1753.11   BIC=1762.02
## 
## Training set error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE        ACF1
## Training set 5.114798 75.44083 52.4759 0.1514047 2.913371 0.6207184 -0.00585451

Significancia 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.458971   0.169486 -2.7080 0.006769 **
## ma1 -0.064159   0.200234 -0.3204 0.748650   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretació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(1,1,1)
## Q* = 86.429, df = 22, p-value = 1.378e-09
## 
## Model df: 2.   Total lags used: 24

Paso 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("Variable1")

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:  86.45698

cat("RMSE Manual: ", rmse_manual, "\n")

## RMSE Manual:  100.7103

MAE (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  2155.514     2096.072
## 2  2024.833  2061.612     2088.980
## 3  2024.917  2085.775     2092.235
## 4  2025.000  2110.975     2090.741
## 5  2025.083  2018.777     2091.427
## 6  2025.167  2174.954     2091.112
## 7  2025.250  2035.961     2091.256
## 8  2025.333  2192.316     2091.190
## 9  2025.417  2118.460     2091.220
## 10 2025.500  2249.288     2091.206
## 11 2025.583  2234.768     2091.213
## 12 2025.667  2207.604     2091.210
## 13 2025.750  2266.162     2091.211
## 14 2025.833  2213.671     2091.211
## 15 2025.917  2218.968     2091.211

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   2096.072
## 2  2024.833   2088.980
## 3  2024.917   2092.235
## 4  2025.000   2090.741
## 5  2025.083   2091.427
## 6  2025.167   2091.112
## 7  2025.250   2091.256
## 8  2025.333   2091.190
## 9  2025.417   2091.220
## 10 2025.500   2091.206
## 11 2025.583   2091.213
## 12 2025.667   2091.210
## 13 2025.750   2091.211
## 14 2025.833   2091.211
## 15 2025.917   2091.211
## 16 2026.000   2091.211

# 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: 2026 = 2091.21076506672"

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: 2096.07184530671"

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,2). 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(4,1,2) with drift 
## 
## Coefficients:
##           ar1     ar2     ar3      ar4      ma1      ma2   drift
##       -0.3909  0.4191  0.0490  -0.2738  -0.1767  -0.4442  3.2478
## s.e.   0.1641  0.1701  0.1216   0.0848   0.1617   0.1337  1.8703
## 
## sigma^2 = 5313:  log likelihood = -864.45
## AIC=1744.89   AICc=1745.9   BIC=1769.09
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.3392507 70.96025 48.58044 -0.1241992 2.710333 0.5746403
##                    ACF1
## Training set 0.01708419

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.390891   0.164082 -2.3823 0.0172050 *  
## ar2    0.419085   0.170053  2.4644 0.0137227 *  
## ar3    0.049049   0.121572  0.4035 0.6866129    
## ar4   -0.273786   0.084835 -3.2273 0.0012497 ** 
## ma1   -0.176674   0.161711 -1.0925 0.2745998    
## ma2   -0.444185   0.133722 -3.3217 0.0008947 ***
## drift  3.247807   1.870313  1.7365 0.0824747 .  
## ---
## 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, 1, 2))  # Especificamos directamente (p=4, d=1, q=2)  

# Mostrar resumen del modelo ajustado
summary(darima_auto)

## Series: train_ts 
## ARIMA(4,1,2) 
## 
## Coefficients:
##           ar1     ar2     ar3      ar4      ma1      ma2
##       -0.4358  0.3604  0.0263  -0.2763  -0.1113  -0.3975
## s.e.   0.1754  0.1850  0.1216   0.0828   0.1747   0.1354
## 
## sigma^2 = 5369:  log likelihood = -865.71
## AIC=1745.43   AICc=1746.2   BIC=1766.59
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 8.781554 71.57468 49.60771 0.342712 2.758073 0.5867916
##                      ACF1
## Training set 0.0008244461

# 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,1,2)
## Q* = 51.552, df = 18, p-value = 4.383e-05
## 
## Model df: 6.   Total lags used: 24

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("variable1")

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

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

# 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  2155.514     2126.319
## 2  2024.833  2061.612     2088.713
## 3  2024.917  2085.775     2126.241
## 4  2025.000  2110.975     2110.058
## 5  2025.083  2018.777     2121.641
## 6  2025.167  2174.954     2126.354
## 7  2025.250  2035.961     2122.184
## 8  2025.333  2192.316     2134.674
## 9  2025.417  2118.460     2128.990
## 10 2025.500  2249.288     2138.837
## 11 2025.583  2234.768     2138.247
## 12 2025.667  2207.604     2142.792
## 13 2025.750  2266.162     2146.693
## 14 2025.833  2213.671     2148.234
## 15 2025.917  2218.968     2153.537

Ahora 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   2126.319
## 2  2024.833   2088.713
## 3  2024.917   2126.241
## 4  2025.000   2110.058
## 5  2025.083   2121.641
## 6  2025.167   2126.354
## 7  2025.250   2122.184
## 8  2025.333   2134.674
## 9  2025.417   2128.990
## 10 2025.500   2138.837
## 11 2025.583   2138.247
## 12 2025.667   2142.792
## 13 2025.750   2146.693
## 14 2025.833   2148.234
## 15 2025.917   2153.537
## 16 2026.000   2154.943

# 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: 2026 = 2154.94331291707"

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: 2126.31875155613"

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,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.4755  0.4279
## s.e.   0.0859  0.0729
## 
## sigma^2 = 4787:  log likelihood = -860.01
## AIC=1726.02   AICc=1726.18   BIC=1735.09

A 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(0, 1, 1),  # (p,d,q) -> (0,1,1)
                seasonal = list(order = c(1, 0, 0),  # (P,D,Q) -> (1,0,0)
                                period = 12))  # Periodicidad estacional de 12 meses

# Mostrar resumen del modelo ajustado
summary(darima)

## Series: train_ts 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.4755  0.4279
## s.e.   0.0859  0.0729
## 
## sigma^2 = 4787:  log likelihood = -860.01
## AIC=1726.02   AICc=1726.18   BIC=1735.09
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 3.099485 68.50635 45.88936 0.06106946 2.564334 0.5428085
##                    ACF1
## Training set 0.02747744

En el correlograma de residuos siguiente se observa que, mejora la correlación de los residuos frente a lso dos modelos anteriores. Sin embargo, al comparar los valores reales VS pronosticados se determina una poca coincidencia. Sigue funcionando mejor el modelo automatico (4,1,2)

# Diagnóstico del modelo (los residuos deben ser ruido blanco)
checkresiduals(darima)  # Verificar si los residuos son aleatorios y no presentan patrones

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(1,0,0)[12]
## Q* = 29.227, df = 22, p-value = 0.1385
## 
## Model df: 2.   Total lags used: 24

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

# 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  2155.514     2098.801
## 2  2024.833  2061.612     2067.649
## 3  2024.917  2085.775     2057.956
## 4  2025.000  2110.975     2064.576
## 5  2025.083  2018.777     2032.200
## 6  2025.167  2174.954     2068.418
## 7  2025.250  2035.961     2051.460
## 8  2025.333  2192.316     2066.896
## 9  2025.417  2118.460     2020.665
## 10 2025.500  2249.288     2066.000
## 11 2025.583  2234.768     2064.054
## 12 2025.667  2207.604     2051.366
## 13 2025.750  2266.162     2059.145
## 14 2025.833  2213.671     2045.816
## 15 2025.917  2218.968     2041.668

Ahora 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   2098.801
## 2  2024.833   2067.649
## 3  2024.917   2057.956
## 4  2025.000   2064.576
## 5  2025.083   2032.200
## 6  2025.167   2068.418
## 7  2025.250   2051.460
## 8  2025.333   2066.896
## 9  2025.417   2020.665
## 10 2025.500   2066.000
## 11 2025.583   2064.054
## 12 2025.667   2051.366
## 13 2025.750   2059.145
## 14 2025.833   2045.816
## 15 2025.917   2041.668
## 16 2026.000   2044.501

# Extraer solo el valor del trimestre adicional (último de la tabla)
next_month <- tail(next_month_forecast, 1)
print(paste("Pronóstico para enero 2025:", next_month$Tiempo, "=", next_month$Pronostico))

## [1] "Pronóstico para enero 2025: 2026 = 2044.5009660821"

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,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.4755  0.4279
## s.e.   0.0859  0.0729
## 
## sigma^2 = 4787:  log likelihood = -860.01
## AIC=1726.02   AICc=1726.18   BIC=1735.09

# 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: 2098.80086926317"

Conclusión:

El modelo automático ARIMA(4,1,2) fue el que mejor desempeño mostró en la comparación entre los datos reales y los pronosticados dentro del periodo de prueba (oct.nov.dic2024). Destacó por su mayor precisión en la captura de los puntos de quiebre, lo que lo hace el más confiable.

No obstante, al analizar los residuos de los modelos, se identifican posibles áreas de mejora para robustecer los pronósticos en los tres casos evaluados. Algunas estrategias podrían incluir la aplicación de una transformación logarítmica o trabajar desde el inicio con la serie ajustada por estacionalidad.

By Julieth Cerón