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("C:/Users/julie/OneDrive - PUJ Cali/Posgrado Finanzas/Maestria Finanzas/Módulo II/Nueva carpeta/Base Caso2.xlsx", 
    col_types = c("date", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))
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
##           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

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.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.86509

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
##            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 1856996.1

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

print(length(tasa_crecimiento_var1))

## [1] 144

print(length(tasa_tendencia_var1))

## [1] 144

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

library(ggplot2)
library(plotly)

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var1)

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

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

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

# Verificar longitudes
print(length(fechas_corregidas_var2))

## [1] 144

print(length(tasa_crecimiento_var2))

## [1] 144

print(length(tasa_tendencia_var2))

## [1] 144

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var2)

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

#Cálculo de la tasa de crecimiento anual correctamente alineada
tasa_crecimiento_var3 <- (variable3_ts[(13:length(variable3_ts))] / variable3_ts[1:(length(variable3_ts) - 12)] - 1) * 100
tasa_tendencia_var3 <- (tendencia_var3[(13:length(tendencia_var3))] / tendencia_var3[1:(length(tendencia_var3) - 12)] - 1) * 100

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

# Verificar longitudes
print(length(fechas_corregidas_var3))

## [1] 144

print(length(tasa_crecimiento_var3))

## [1] 144

print(length(tasa_tendencia_var3))

## [1] 144

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

# Convertir a gráfico interactivo
ggplotly(grafico_crecimiento_var3)

Analizar la tasa de crecimiento anual ayuda a detectar cambios en el entorno económico que afectan el sector. Se pueden prever crisis o períodos de auge y prepararse para ellos.

Modelo ARIMA

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

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

## RMSE Manual:  37.96552

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

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

# 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 = 2090.74082063412"

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

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

# 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 = 2110.05769478518"

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

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

# 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 = 2064.57613893857"

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