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

library(readxl)
data_col <- read_excel("/Users/sofiaartunduaga/Desktop/Econometría/Históricos BRNT.xlsx")
print(data_col)
## # A tibble: 120 × 7
##    Fecha               Último Apertura Máximo Mínimo Vol.   `% var.`
##    <dttm>               <dbl>    <dbl>  <dbl>  <dbl> <chr>     <dbl>
##  1 2015-04-01 00:00:00   38.6     33.2   39.0   32.8 21,85K   0.171 
##  2 2015-05-01 00:00:00   37.2     38.8   40.5   35.4 75,06K  -0.0354
##  3 2015-06-01 00:00:00   35.9     37.3   38.2   35.4 4,45K   -0.0359
##  4 2015-07-01 00:00:00   30.1     35.8   36.2   30.0 3,47K   -0.161 
##  5 2015-08-01 00:00:00   28.0     29.5   29.5   24.4 35,29K  -0.0692
##  6 2015-09-01 00:00:00   26.7     29.8   30.0   26.2 1,51K   -0.0482
##  7 2015-10-01 00:00:00   26.8     27.2   29.5   25.3 0,53K    0.0052
##  8 2015-11-01 00:00:00   24.2     26.5   27.7   23.4 0,15K   -0.0981
##  9 2015-12-01 00:00:00   19.1     23.6   23.7   19.0 3,80K   -0.210 
## 10 2016-01-01 00:00:00   17.8     19.9   20.6   14.2 52,91K  -0.0715
## # ℹ 110 more rows
# Convertir/declarar variable 1=ENER en serie de tiempo mensual
variable1_ts <- ts(data_col$Último, start = c(2015, 1), frequency = 12)
variable1_ts
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2015 38.60 37.24 35.90 30.11 28.03 26.68 26.82 24.19 19.12 17.75 18.00 19.60
## 2016 22.51 23.86 23.45 19.92 21.56 22.50 21.82 22.10 24.23 23.80 23.40 22.22
## 2017 21.64 20.82 19.77 21.19 21.23 23.13 24.77 25.67 27.23 28.27 27.19 28.84
## 2018 31.15 32.84 33.53 31.89 33.19 35.46 32.82 25.27 23.44 26.58 28.36 28.86
## 2019 31.03 27.63 28.87 28.64 26.17 27.15 27.19 28.33 30.98 27.40 23.86 13.53
## 2020 12.12 15.52 17.85 18.49 19.43 17.38 15.55 19.23 20.68 22.41 26.50 26.44
## 2021 27.54 28.77 31.00 31.60 30.58 33.59 36.03 31.02 34.39 39.33 43.80 50.56
## 2022 51.97 56.64 54.44 53.94 50.35 45.70 49.96 48.24 46.46 47.43 46.83 44.79
## 2023 45.44 42.02 43.44 49.40 50.20 53.99 51.72 48.26 46.96 49.47 50.63 53.63
## 2024 54.35 51.81 54.33 52.34 50.69 48.01 48.73 49.03 50.46 52.18 50.58 49.22

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)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# 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("2015-04-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)

Gráfica 1

Este gráfico muestra la evolución del precio del WisdomTree Brent Crude Oil ETC a lo largo de la década, revelando importantes altibajos. Se identifican tres períodos clave:

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

Gráfica 2

Esta gráfica muestra la descomposición temporal del WisdomTree Brent Crude Oil ETC, separando la serie original en sus distintos componentes: tendencia, estacionalidad y residuo. La serie original representa el comportamiento real de los precios del petróleo entre 2015 y 2025, reflejando períodos de caída, recuperación y estabilidad. Se observa una fuerte caída entre 2015 y 2016, una recuperación progresiva hasta 2019 y una abrupta disminución en 2020 debido a la pandemia de COVID-19. Posteriormente, el precio muestra un crecimiento hasta 2022 y luego una fase de estabilización hasta 2025, lo que indica que el mercado ha encontrado un punto de equilibrio.

La tendencia revela la dirección general del mercado, mostrando una caída sostenida hasta 2020, seguida de una recuperación acelerada y una posterior estabilización. La estacionalidad refleja fluctuaciones predecibles que pueden estar relacionadas con la demanda de petróleo en diferentes épocas del año, como el invierno o el verano. Aunque la demanda de petróleo suele experimentar variaciones estacionales, como el aumento en invierno debido al mayor uso de calefacción o el incremento en verano por el mayor consumo de combustible en viajes, en este caso, la influencia de estos patrones es mínima en comparación con los eventos globales que han determinado la evolución del precio. Factores como crisis económicas, conflictos geopolíticos y cambios en la producción han tenido un impacto mucho más significativo en las fluctuaciones del mercado, opacando cualquier efecto estacional predecible. Por otro lado, el residuo recoge variaciones impredecibles no explicadas por la tendencia o la estacionalidad, destacando picos abruptos en momentos clave como la crisis del COVID-19 en 2020 o el conflicto entre Rusia y Ucrania en 2022.

Se crea la variable1 ajustada por estacionalidad

# Extraer los componentes de la descomposición
variable1_sa <- variable1_ts - stl_decomp_var1$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("2015-04-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 = "red", size = 0.5, linetype = "solid", name = "Serie Original") +
  geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", size = 0.6, linetype = "solid", name = "Serie Ajustada") +
  ggtitle("Variable 1:Serie Original vs Serie Ajustada por Estacionalidad") +
  xlab("Tiempo") +
  ylab("Unidad de medida variable 1") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor visualización
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning in geom_line(aes(x = fechas_var1, y = variable1_ts), color = "red", :
## Ignoring unknown parameters: `name`
## Warning in geom_line(aes(x = fechas_var1, y = variable1_sa), color = "black", :
## Ignoring unknown parameters: `name`
# Convertir a gráfico interactivo
ggplotly(grafico_ajustada_var1)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Gráfico 3

Este gráfico compara la serie original con una versión ajustada que elimina los efectos estacionales, permitiendo observar con mayor claridad la dirección del mercado. Se puede notar que ambas curvas son muy similares, lo que confirma que los cambios de precio se deben más a eventos globales y tendencias de largo plazo que a fluctuaciones estacionales.

La importancia de esta comparación radica en que ayuda a los inversionistas a distinguir entre movimientos temporales y cambios estructurales en el precio del petróleo. Al eliminar la estacionalidad, es más fácil prever la dirección futura del mercado sin dejarse influenciar por cambios pasajeros.

Primero se debe obtener la tendencia de cada variable y luego graficarla

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("2015-04-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" = "blue", "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)

Gráfico 4

Comparar la serie original con su tendencia nos ayuda a distinguir entre la volatilidad a corto plazo y la dirección general del mercado.

En este caso, la tendencia confirma que entre 2015 y 2020 hubo una fase bajista, con una disminución progresiva del precio. Luego, a partir del 2020, se observa un crecimiento fuerte hasta 2022, seguido de una estabilización. Esto indica que, aunque el mercado experimenta muchas fluctuaciones, en el largo plazo hay períodos bien definidos de caída, recuperación y estabilidad.

Para los inversionistas, la tendencia es una herramienta clave, ya que permite analizar el mercado sin dejarse llevar por cambios repentinos que pueden ser engañosos.

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 abril 2015
fechas_corregidas_var1 <- seq(from = as.Date("2015-04-01"), by = "month", length.out = length(tasa_crecimiento_var1))

# Verificar longitudes
print(length(fechas_corregidas_var1))
## [1] 108
print(length(tasa_crecimiento_var1))
## [1] 108
print(length(tasa_tendencia_var1))
## [1] 108
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 = "black", size = 0.7) +
  geom_line(aes(x = fechas_corregidas_var1, y = tasa_tendencia_var1), color = "red", 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)

Gráfico 5

El último gráfico nos muestra cómo han cambiado los precios en términos porcentuales año tras año. Aquí se hace evidente que el mercado del petróleo sigue un comportamiento cíclico: después de cada caída importante, hay un período de recuperación.

Por ejemplo, tras la fuerte caída de 2015-2016, hubo una fase de crecimiento entre 2016 y 2018. Luego, en 2018-2020, el precio volvió a caer, alcanzando su punto más bajo en 2020. Sin embargo, después de este mínimo, se dio un fuerte repunte entre 2020 y 2022, con tasas de crecimiento muy elevadas. Finalmente, entre 2022 y 2025, el crecimiento se ha estabilizado, con menos fluctuaciones extremas.

Este comportamiento se explica por la dinámica de oferta y demanda. Cuando los precios caen demasiado, las empresas petroleras reducen la producción porque deja de ser rentable. Con menos oferta en el mercado, los precios comienzan a subir nuevamente. A esto se suman factores como la recuperación económica y decisiones estratégicas de la OPEP, que regulan la producción para evitar caídas prolongadas.

Este patrón cíclico es clave para quienes invierten en petróleo, por ejemplo, ya que les permite anticipar momentos de compra y venta según la fase del ciclo en la que se encuentre el mercado.

Conclusión

El comportamiento del WisdomTree Brent Crude Oil ETC entre 2015 y 2025 refleja la interacción de múltiples factores que han moldeado la evolución del mercado del petróleo. Desde crisis económicas hasta tensiones geopolíticas, cada evento ha dejado una huella significativa en los precios. La estabilización reciente sugiere que el mercado ha encontrado un nuevo equilibrio, aunque sigue siendo vulnerable a futuros cambios en la oferta y la demanda global. A partir de 2025, se espera una relativa estabilidad en torno a las 45-50 unidades, aunque la creciente transición hacia energías renovables y las regulaciones ambientales podrían generar una presión bajista a largo plazo. Sin embargo, eventos geopolíticos o disrupciones en la oferta pueden provocar aumentos inesperados. En este contexto, las empresas de transporte y manufactura, que dependen del petróleo, se benefician de precios más bajos al reducir sus costos operativos, al igual que las economías importadoras de crudo, como las europeas y asiáticas, que pueden mejorar su balanza comercial. En contraste, los países altamente dependientes de la exportación de petróleo, como Rusia y los productores de Medio Oriente, podrían enfrentar dificultades si la demanda disminuye o los precios se mantienen bajos. Asimismo, las petroleras tradicionales deben adaptarse a una posible reducción de márgenes de ganancia ante la aceleración de la transición energética. Esto resalta la importancia de comprender tanto los movimientos cíclicos como los estructurales para anticipar tendencias futuras y tomar decisiones informadas en el mercado energético.

# 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
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# 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 = -2.3259, Lag order = 4, p-value = 0.4413
## alternative hypothesis: stationary
#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) 
library(plotly)
 # 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
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Interpretación

  # 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
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
 # 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)

library(tseries)
# Segunda prueba de estacionariedad sobre la serie diferenciada en niveles
  adf_test_diff <- adf.test(train_diff)
## Warning in adf.test(train_diff): p-value smaller than printed p-value
  print(adf_test_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff
## Dickey-Fuller = -4.8222, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
# Segunda prueba de estacionariedad sobre la serie diferenciada en logaritmo
  adf_test_diff_log <- adf.test(train_diff_log)
## Warning in adf.test(train_diff_log): p-value smaller than printed p-value
  print(adf_test_diff_log)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff_log
## Dickey-Fuller = -4.9697, Lag order = 4, 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(ggplot2)
library(plotly)
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)

Interpretación Este gráfico muestra la función de autocorrelación de la serie temporal antes de aplicar el modelo ARIMA. Sirve para identificar la estructura de los datos y seleccionar el valor de q en el modelo ARIMA(p,d,q).

Se observa un pico significativo en el lag 1, lo que sugiere que la serie presenta una memoria corta. Esto significa que el valor actual de la serie está fuertemente correlacionado con su valor inmediato anterior. También hay otros valores menores, pero la caída rápida de la autocorrelación indica que un modelo MA(1) (media móvil de orden 1) podría ser adecuado para capturar esta dependencia.

Si los valores de la autocorrelación disminuyeran gradualmente en vez de caer abruptamente, sería una señal de que un modelo AR (autorregresivo) sería más apropiado. En este caso, el comportamiento observado sugiere que la inclusión de un término de media móvil q=1 podría mejorar el ajuste del modelo.

ggplotly(pacf_plot)

Interpretación La función de autocorrelación parcial (PACF) es utilizada para determinar el orden p del modelo ARIMA, es decir, cuántos valores pasados de la serie influyen directamente en el presente.

En este gráfico se observa un pico significativo en el lag 1, lo que indica que la serie tiene un componente autorregresivo fuerte en ese punto. También hay algunos valores menores en otros lags, pero el hecho de que la correlación caiga después del primer lag sugiere que un modelo AR(1) (autorregresivo de orden 1) podría ser suficiente.

El análisis conjunto de ACF y PACF indica que el modelo ARIMA(1,1,1) es una opción razonable, ya que se identificaron tanto un componente autorregresivo (p=1) como de media móvil (q=1). Sin embargo, la presencia de ciertos valores significativos en ambos gráficos sugiere que se podrían probar otros modelos, como ARIMA(1,1,2) o ARIMA(2,1,1), para mejorar la captura de la estructura de la serie.

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.5113  0.7462
## s.e.   0.2634  0.2136
## 
## sigma^2 = 6.804:  log likelihood = -274.88
## AIC=555.76   AICc=555.98   BIC=564.02
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.09191719 2.574777 1.988331 -0.1908026 6.774285 0.2374005
##                    ACF1
## Training set 0.05101358

Significancia de coefientes

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(manual_arima_model)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.51128    0.26344 -1.9408 0.0522847 .  
## ma1  0.74621    0.21362  3.4932 0.0004773 ***
## ---
## 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* = 22.969, df = 21, p-value = 0.3456
## 
## Model df: 2.   Total lags used: 23

Interpretación Serie de Residuos a lo largo del tiempo: Muestra la evolución de los residuos del modelo. Se espera que los residuos sean aleatorios, con media cero y sin una tendencia clara. Si hay patrones visibles, puede indicar que el modelo no está capturando completamente la estructura de los datos. Función de Autocorrelación - ACF: Muestra la autocorrelación de los residuos en diferentes retardos (lags). Las líneas punteadas azules representan los intervalos de confianza. Si los valores caen dentro de estas líneas, no hay autocorrelación significativa. En este caso, hay algunos picos que podrían indicar que los residuos no son completamente aleatorios. Histograma de los residuos con ajuste de densidad: Evalúa si los residuos siguen una distribución normal. Si el histograma es simétrico y se ajusta bien a la curva normal (línea naranja), los residuos pueden considerarse normalmente distribuidos. Aquí parece haber una forma cercana a la normal, pero con una ligera asimetría.

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)
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Interpretar En la gráfica vemos dos líneas: una representa los valores reales de la serie observada y la otra muestra el pronóstico manual realizado. Lo primero que llama la atención es que el pronóstico manual no sigue del todo el comportamiento de los valores reales. Mientras la serie observada tiene una tendencia a la baja en varios periodos, el pronóstico manual es más estable y hasta muestra ligeras subidas en algunos puntos. Esto indica que el pronóstico manual no está capturando correctamente la tendencia real del fenómeno estudiado, lo que podría llevar a decisiones erróneas basadas en una expectativa incorrecta del futuro.

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:  1.006399
cat("RMSE Manual: ", rmse_manual, "\n")
## RMSE Manual:  1.19042

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)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# 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     52.18     51.21027
## 2 2024.833     50.58     50.82667
## 3 2024.917     49.22     51.02280

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   51.21027
## 2 2024.833   50.82667
## 3 2024.917   51.02280
## 4 2025.000   50.92252
# 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 = 50.9225201593705"

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

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

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(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.2646
## s.e.  0.0923
## 
## sigma^2 = 6.767:  log likelihood = -275.03
## AIC=554.06   AICc=554.17   BIC=559.57
## 
## Training set error measures:
##                      ME     RMSE     MAE        MPE     MAPE      MASE
## Training set 0.08331647 2.579044 1.98406 -0.1500843 6.752849 0.2368906
##                      ACF1
## Training set 0.0009360508
library(lmtest)

# Evaluar la significancia estadística de los coeficientes del modelo ARIMA
coeftest(auto_arima_model_no_seasonal)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)   
## ma1 0.264618   0.092279  2.8676 0.004136 **
## ---
## 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.2871  0.1275  -0.0177  -0.1852  0.5697  -0.0343
## s.e.   0.4128  0.4391   0.1538   0.1094  0.4213   0.4566
## 
## sigma^2 = 6.765:  log likelihood = -272.54
## AIC=559.07   AICc=560.11   BIC=578.35
## 
## Training set error measures:
##                      ME    RMSE      MAE        MPE     MAPE      MASE
## Training set 0.09498328 2.52197 1.980032 -0.1482074 6.766732 0.2364096
##                      ACF1
## Training set 0.0001476243
# 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* = 17.771, df = 17, p-value = 0.4034
## 
## Model df: 6.   Total lags used: 23

Interpretar En este caso, se probó un modelo ARIMA(4,1,2), lo que significa que está diseñado para captar patrones de tendencia y fluctuaciones en los datos.

Para evaluar el desempeño, se analizaron varias métricas de error. El error medio (ME) es cercano a cero, lo que sugiere que el modelo no tiene un sesgo fuerte (es decir, no tiende a sobrestimar ni subestimar de manera constante). El RMSE y el MAE nos indican el tamaño promedio del error en las predicciones, y en este caso, aunque no son valores perfectos, están dentro de un rango aceptable.

Otro análisis importante es el de los residuos, que básicamente son las diferencias entre los valores reales y los valores que el modelo predijo. Si los residuos se comportan de manera aleatoria y no tienen patrones, significa que el modelo hizo un buen trabajo capturando la estructura de los datos. Aquí, los residuos se distribuyen alrededor de cero sin mostrar tendencias claras, lo que es una señal positiva. Además, al analizar la autocorrelación de los residuos (ACF), no se encontraron correlaciones significativas, lo que refuerza la idea de que el modelo está bien ajustado.

Pronóstico modelo ARIMA automático (4,1,2)

library(plotly)

# 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
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
# 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     52.18     50.84271
## 2 2024.833     50.58     50.84271
## 3 2024.917     49.22     50.84271

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   50.84271
## 2 2024.833   50.84271
## 3 2024.917   50.84271
## 4 2025.000   50.84271
# 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 = 50.8427111965766"
# 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: 50.8427111965766"

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)(0,0,1)[12] 
## 
## Coefficients:
##          ma1     sma1
##       0.2516  -0.1858
## s.e.  0.0939   0.0988
## 
## sigma^2 = 6.61:  log likelihood = -273.37
## AIC=552.73   AICc=552.95   BIC=560.99

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.2533  -0.1601
## s.e.  0.0932   0.0941
## 
## sigma^2 = 6.644:  log likelihood = -273.61
## AIC=553.22   AICc=553.43   BIC=561.48
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.09398089 2.544281 1.934377 -0.1741098 6.597164 0.2309586
##                      ACF1
## Training set -0.002082818

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* = 16.69, df = 21, p-value = 0.7297
## 
## Model df: 2.   Total lags used: 23

Interpretación

El modelo SARIMA es una versión mejorada del ARIMA que incorpora componentes estacionales, es decir, patrones que se repiten cada cierto tiempo. En este caso, se utilizó un modelo SARIMA con periodicidad de 12 unidades, lo que sugiere que se esperaba capturar algún tipo de repetición anual en los datos.

Las métricas de error de este modelo son muy similares a las del ARIMA, aunque con una ligera ventaja en términos de MAPE (error porcentual absoluto medio). En cuanto al análisis de residuos, este modelo mostró una menor correlación en los errores en comparación con el ARIMA, lo cual es un punto a favor, ya que indica que está eliminando más estructura de los datos y dejando residuos más cercanos a un comportamiento aleatorio.

Sin embargo, la diferencia en desempeño entre ambos modelos no es muy grande, lo que significa que el modelo ARIMA sigue siendo una opción competitiva.

Pronóstico con el modelo SARIMA

library(ggplot2)
# Generar pronóstico para el conjunto de prueba
forecast_arima <- forecast(darima, h = length(test_ts))  # Predecir los valores futuros

# Crear dataframe para gráfico interactivo del pronóstico
forecast_data <- data.frame(Tiempo = time(forecast_arima$mean), 
                            Pronostico = as.numeric(forecast_arima$mean),
                            Observado = as.numeric(test_ts))

# Graficar pronóstico junto con los valores observados reales
p4 <- ggplot(forecast_data, aes(x = Tiempo)) +
  geom_line(aes(y = Pronostico, color = "Pronóstico")) +
  geom_line(aes(y = Observado, color = "Observado")) +
  ggtitle("Pronóstico vs Observado") +
  xlab("Tiempo") + ylab("Unidad Variable 1")

ggplotly(p4)  # Convertir el gráfico en interactivo
## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.

Interpretación

Por último, se comparó el pronóstico generado automáticamente por los modelos con la serie observada. En este caso, el modelo automático muestra una mejor aproximación a la realidad que el pronóstico manual. Aunque no es perfecto, se puede ver que sigue más de cerca las tendencias y fluctuaciones de la serie original.

Esto nos dice que utilizar modelos estadísticos para hacer predicciones puede ser mucho más efectivo que hacer pronósticos basados en intuición o métodos simples. Sin embargo, siempre hay margen de mejora, y dependiendo de la aplicación específica, se podrían probar modelos más avanzados para reducir aún más el error.

# 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     52.18     50.25844
## 2 2024.833     50.58     50.03079
## 3 2024.917     49.22     49.50277

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   50.25844
## 2 2024.833   50.03079
## 3 2024.917   49.50277
## 4 2025.000   49.33641
# 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 = 49.336405720423"

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)(0,0,1)[12] 
## 
## Coefficients:
##          ma1     sma1
##       0.2516  -0.1858
## s.e.  0.0939   0.0988
## 
## sigma^2 = 6.61:  log likelihood = -273.37
## AIC=552.73   AICc=552.95   BIC=560.99
# 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: 50.2584380996599"

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.