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

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.

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"