Análisis y Pronóstico de Precios del Petróleo Mexicano Utilizando Modelos ARIMA

# Cargar paquetes necesarios
packages <- c("readxl", "ggplot2", "forecast", "dplyr", "tseries", "plotly")
lapply(packages, function(pkg) if (!require(pkg, character.only = TRUE)) install.packages(pkg))

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL

# Cargar las librerías
library(readxl)
library(ggplot2)
library(forecast)
library(dplyr)
library(tseries)
library(plotly)
# Leer el archivo Excel
archivo <- "C:/Users/MARIN/Documents/Precios_2.xlsx"
datos <- read_excel(archivo)


# Verificar las primeras filas para asegurarse de que se cargó correctamente
head(datos)

## # A tibble: 6 × 2
##   Fecha               Precios           
##   <dttm>              <chr>             
## 1 1996-01-03 00:00:00 17.399999999999999
## 2 1996-01-04 00:00:00 17.41             
## 3 1996-01-05 00:00:00 17.7              
## 4 1996-01-08 00:00:00 17.54             
## 5 1996-01-09 00:00:00 17.41             
## 6 1996-01-10 00:00:00 17.21

# Reemplazar "N/E" por NA en la columna de precios
# Asume que la columna de precios se llama 'Precios'
datos <- datos %>%
  mutate(Precios = 
           na_if(Precios, "N/E")) %>%
  mutate(Precios =
           as.numeric(Precios))  # Convertir a numérico

# Verificar las primeras filas para asegurarse de que se limpió correctamente
head(datos)

## # A tibble: 6 × 2
##   Fecha               Precios
##   <dttm>                <dbl>
## 1 1996-01-03 00:00:00    17.4
## 2 1996-01-04 00:00:00    17.4
## 3 1996-01-05 00:00:00    17.7
## 4 1996-01-08 00:00:00    17.5
## 5 1996-01-09 00:00:00    17.4
## 6 1996-01-10 00:00:00    17.2

# Filtrar los datos desde 2022
datos_filtrados <- subset(datos, Fecha >= as.Date("2018-12-03"))[, c("Fecha", "Precios")]

# Crear el gráfico de líneas
ggplot(datos_filtrados, aes(x = Fecha, y = Precios)) +
  geom_line(color = "blue") +
  labs(title = "Precio del Petróleo (USD) desde 2018", x = "Fecha", y = "Precio en USD") +
  theme_minimal()

#base de datos descargada de url:https://www.banxico.org.mx/apps/gc/precios-spot-del-petroleo-gra.html

# Filtrar los datos desde 2022
datos_filtrados <- subset(datos, Fecha >= as.Date("2018-12-03"))[, c("Fecha", "Precios")]

# Crear el gráfico de líneas con una línea de tendencia
ggplot(datos_filtrados, aes(x = Fecha, y = Precios)) +
  geom_line(color = "blue") +  # Línea de los precios
  geom_smooth(method = "lm", se = FALSE, color = "red") +  # Línea de tendencia
  labs(title = "Precio del Petróleo (USD) desde 2018 con Línea de Tendencia", x = "Fecha", y = "Precio en USD") +
  theme_minimal()

#Con el grafico anterior podemos observar que no es estacionario y presenta tendencia 
#Por lo anterior procederemos a aplicarle una primera diferencia con el fin que nuestra serie se comporte más estacionaria 

# Filtrar los datos desde 2022
datos_filtrados <- subset(datos, Fecha >= as.Date("2018-12-03"))[, c("Fecha", "Precios")]

# Calcular la primera diferencia de la columna de precios
datos_filtrados$Diferencia_Precio <- c(NA, diff(datos_filtrados$Precios))

# Crear el gráfico de la primera diferencia
ggplot(datos_filtrados, aes(x = Fecha, y = Diferencia_Precio)) +
  geom_line(color = "grey") +
  labs(title = "Primera Diferencia del Precio del Petróleo (USD)", x = "Fecha", y = "Diferencia del Precio") +
  theme_minimal()

#Se observa que la serie fluctúa alrededor de una media estable y finita salvo lo que parece ser el 2011 

# Asegurarse de que la columna de precios es numérica
datos$Precios <- as.numeric(datos$Precios)

# Eliminar valores NA
precios_limpios <- na.omit(datos$Precios)

# Verificar que la serie no tenga valores constantes
if (length(unique(precios_limpios)) > 1) {
  # Realizar la prueba de raíz unitaria (Dickey-Fuller) solo si la serie no es constante
  adf_resultado <- adf.test(precios_limpios, alternative = "stationary")
  
  # Mostrar el resultado de la prueba
  print(adf_resultado)
} else {
  print("La serie de precios es constante o tiene pocos valores únicos, no se puede realizar la prueba de ADF.")
}

## 
##  Augmented Dickey-Fuller Test
## 
## data:  precios_limpios
## Dickey-Fuller = -2.5659, Lag order = 19, p-value = 0.3387
## alternative hypothesis: stationary

#dado el valor p no podemos rechazar la hipotesis nula lo que sugiere que la serie no es estacionaria  y procedermos a sacar la 2a diferencia

# Filtrar los datos desde 2022
datos_filtrados <- subset(datos, Fecha >= as.Date("2018-12-03"))[, c("Fecha", "Precios")]

# Calcular la segunda diferencia de la columna de precios
datos_filtrados$Segunda_Diferencia_Precio <- c(NA, NA, diff(diff(datos_filtrados$Precios)))

# Crear el gráfico de la segunda diferencia
ggplot(datos_filtrados, aes(x = Fecha, y = Segunda_Diferencia_Precio)) +
  geom_line(color = "orange") +
  labs(title = "Segunda Diferencia del Precio del Petróleo (USD)", x = "Fecha", y = "Segunda Diferencia del Precio") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Realizar la prueba de raíz unitaria (Dickey-Fuller) en la segunda diferencia de precios
adf_segunda_diferencia <- adf.test(na.omit(datos_filtrados$Segunda_Diferencia_Precio), alternative = "stationary")

# Mostrar el resultado de la prueba
print(adf_segunda_diferencia)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  na.omit(datos_filtrados$Segunda_Diferencia_Precio)
## Dickey-Fuller = -21.008, Lag order = 18, p-value = 0.01
## alternative hypothesis: stationary

#Interpretación del resultado de la prueba de Dickey-Fuller aumentada (ADF)

#Resultado de la prueba:
#Dickey-Fuller = -21.008: Este es el valor estadístico de la prueba. Un valor más negativo indica una mayor evidencia contra la hipótesis nula (existencia de una raíz unitaria, o no estacionariedad). En este caso, -21.008 es un valor muy negativo, lo que sugiere que hay fuerte evidencia en contra de la hipótesis nula.

#Lag order = 18: Indica que se utilizaron 18 retardos (lags) en el modelo. Este número se selecciona para capturar la autocorrelación en los datos.

#p-value = 0.01: Este es el valor p de la prueba. Un valor p de 0.01 es bastante bajo, lo que indica que hay suficiente evidencia para rechazar la hipótesis nula de que la serie tiene una raíz unitaria. Esto sugiere que la serie es estacionaria.

#alternative hypothesis: stationary: Esto indica que la hipótesis alternativa es que la serie es estacionaria (es decir, no tiene una raíz unitaria).

#Interpretación:
#Hipótesis nula (H0): La serie tiene una raíz unitaria (no es estacionaria).
#Hipótesis alternativa (H1): La serie es estacionaria (no tiene una raíz unitaria).
#Conclusión:
#Dado que el valor p es 0.01, que es menor que el umbral típico de 0.05, puedes rechazar la hipótesis nula. Esto sugiere que la serie de la segunda diferencia de los precios del petróleo es estacionaria. Esto es una buena noticia.

#Para calcular y visualizar las funciones de autocorrelación (ACF) y autocorrelación parcial (PACF) de la serie de precios o de la segunda diferencia, uso las funciones acf() y pacf() en R. Esto  ayudará a identificar la naturaleza del proceso generador de datos y a determinar los componentes de los modelos ARIMA.


# Calcular y graficar la función de autocorrelación (ACF)
acf_result <- acf(na.omit(datos_filtrados$Segunda_Diferencia_Precio), 
                  main = "Función de Autocorrelación (ACF)", 
                  lag.max = 40,  # Puedes ajustar el número de lags a mostrar
                  plot = TRUE)

# Calcular y graficar la función de autocorrelación parcial (PACF)
pacf_result <- pacf(na.omit(datos_filtrados$Segunda_Diferencia_Precio), 
                    main = "Función de Autocorrelación Parcial (PACF)", 
                    lag.max = 40,  # Puedes ajustar el número de lags a mostrar
                    plot = TRUE)

#con base en los resultados anteriores,pronosticaremos con ARIMA para AR(1,2,3,4,5,6,7,8) Y MA(1,2)
# Estimar el mejor modelo ARIMA automáticamente
modelo_arima_auto <- auto.arima(na.omit(datos_filtrados$Segunda_Diferencia_Precio))

# Mostrar el resumen del modelo
summary(modelo_arima_auto)

## Series: na.omit(datos_filtrados$Segunda_Diferencia_Precio) 
## ARIMA(0,0,2) with zero mean 
## 
## Coefficients:
##           ma1      ma2
##       -0.6667  -0.0452
## s.e.   0.0130   0.0133
## 
## sigma^2 = 2.358:  log likelihood = -11293.83
## AIC=22593.67   AICc=22593.67   BIC=22613.82
## 
## Training set error measures:
##                      ME     RMSE       MAE MPE MAPE      MASE          ACF1
## Training set 0.01049389 1.535444 0.9598601 NaN  Inf 0.4854334 -0.0009369964

# Inicializar una lista para almacenar los modelos
resultados <- list()

# Definir los órdenes AR y MA a probar
ar_orders <- c(1, 2, 3, 4, 5, 6, 7, 8, 11, 22, 25)
ma_orders <- c(1, 2, 3, 5, 6)

# Contador para limitar a 5 combinaciones
contador <- 0

# Ajustar modelos ARIMA hasta 5 combinaciones
for (p in ar_orders) {
  for (q in ma_orders) {
    if (contador < 5) {
      # Intentar ajustar el modelo ARIMA
      modelo_arima <- tryCatch({
        arima(na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(p, 0, q))
      }, error = function(e) {
        return(NULL)  # Si hay un error, devuelve NULL
      })
      
      # Almacenar el modelo solo si fue exitoso
      if (!is.null(modelo_arima)) {
        resultados[[paste("AR", p, "MA", q, sep = "_")]] <- modelo_arima
        contador <- contador + 1  # Aumentar el contador
      }
    } else {
      break  # Salir del bucle si ya se han ajustado 5 modelos
    }
  }
  if (contador >= 5) break  # Salir del bucle si ya se han ajustado 5 modelos
}

# Mostrar un resumen de los modelos ajustados
model_summaries <- lapply(resultados, summary)

# Mostrar un ejemplo de los resultados (opcional)
model_summaries

## $AR_1_MA_1
## 
## Call:
## arima(x = na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(p, 0, 
##     q))
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.0615  -0.7299     0.0030
## s.e.  0.0185   0.0128     0.0057
## 
## sigma^2 estimated as 2.358:  log likelihood = -11293.98,  aic = 22595.96
## 
## Training set error measures:
##                        ME    RMSE       MAE MPE MAPE     MASE         ACF1
## Training set 2.372122e-07 1.53548 0.9596873 NaN  Inf 0.485346 0.0006024604
## 
## $AR_1_MA_2
## 
## Call:
## arima(x = na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(p, 0, 
##     q))
## 
## Coefficients:
##           ar1      ma1      ma2  intercept
##       -0.2680  -0.3981  -0.2329     0.0030
## s.e.   0.1607   0.1591   0.1097     0.0057
## 
## sigma^2 estimated as 2.357:  log likelihood = -11292.71,  aic = 22595.42
## 
## Training set error measures:
##                        ME     RMSE       MAE MPE MAPE      MASE         ACF1
## Training set -1.18367e-05 1.535162 0.9594952 NaN  Inf 0.4852488 -0.001321625
## 
## $AR_1_MA_3
## 
## Call:
## arima(x = na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(p, 0, 
##     q))
## 
## Coefficients:
##           ar1      ma1      ma2     ma3  intercept
##       -0.2083  -0.4588  -0.1952  0.0081     0.0030
## s.e.   0.2256   0.2253   0.1494  0.0166     0.0058
## 
## sigma^2 estimated as 2.357:  log likelihood = -11292.57,  aic = 22597.14
## 
## Training set error measures:
##                        ME     RMSE       MAE MPE MAPE      MASE          ACF1
## Training set 2.058885e-05 1.535127 0.9593955 NaN  Inf 0.4851984 -0.0002562498
## 
## $AR_1_MA_5
## 
## Call:
## arima(x = na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(p, 0, 
##     q))
## 
## Coefficients:
##           ar1      ma1      ma2     ma3     ma4      ma5  intercept
##       -0.2608  -0.4057  -0.2232  0.0258  0.0251  -0.0792     0.0030
## s.e.   0.1165   0.1158   0.0782  0.0152  0.0141   0.0128     0.0053
## 
## sigma^2 estimated as 2.342:  log likelihood = -11273.38,  aic = 22562.76
## 
## Training set error measures:
##                        ME     RMSE       MAE MPE MAPE      MASE          ACF1
## Training set 8.313546e-07 1.530312 0.9601296 NaN  Inf 0.4855697 -0.0008369686
## 
## $AR_1_MA_6
## 
## Call:
## arima(x = na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(p, 0, 
##     q))
## 
## Coefficients:
##          ar1      ma1     ma2     ma3      ma4      ma5     ma6  intercept
##       0.8275  -1.4952  0.5009  0.0805  -0.0190  -0.1004  0.0868     0.0031
## s.e.  0.0716   0.0724  0.0528  0.0240   0.0243   0.0235  0.0131     0.0061
## 
## sigma^2 estimated as 2.338:  log likelihood = -11267.88,  aic = 22553.76
## 
## Training set error measures:
##                        ME     RMSE       MAE MPE MAPE      MASE          ACF1
## Training set -6.80553e-05 1.528931 0.9590009 NaN  Inf 0.4849989 -0.0009513372

#Interpretación de Resultados
#Modelos ARIMA Estimados:

#$AR_1_MA_1: ARIMA(1, 0, 1)
#$AR_1_MA_2: ARIMA(1, 0, 2)
#$AR_1_MA_3: ARIMA(1, 0, 3)
#$AR_1_MA_5: ARIMA(1, 0, 5)
#$AR_1_MA_6: ARIMA(1, 0, 6)
#Coeficientes:

#Cada modelo muestra diferentes coeficientes para el componente autorregresivo (AR) y para los componentes de media móvil (MA). Por ejemplo:
#En el modelo AR_1_MA_6, el coeficiente AR es positivo (0.8275), lo que sugiere que el valor actual de la serie tiene una relación positiva con su valor anterior.
#Los coeficientes MA pueden ser negativos o positivos, lo que indica la relación entre el error en la predicción y los valores pasados.
#Errores Cuadráticos Medios (RMSE):

#Los valores de RMSE son bastante similares entre los modelos, indicando que todos los modelos tienen un error de ajuste similar, pero disminuyen ligeramente con el aumento del orden de MA.
#AR_1_MA_6 tiene el RMSE más bajo (1.528931), lo que sugiere que podría ser el mejor modelo en términos de ajuste.
#Criterios de Información:

#El AIC (Criterio de Información de Akaike) también es un buen criterio para comparar modelos. El modelo AR_1_MA_6 tiene el AIC más bajo (22553.76), lo que es favorable.
#En general, un AIC más bajo indica un mejor ajuste con menos complejidad.
#Log Likelihood:

#Un valor de log likelihood más alto (menos negativo) indica un mejor ajuste al modelo. AR_1_MA_6 tiene el valor de log likelihood más alto (-11267.88), lo que respalda su elección como el mejor modelo.
#Mejor Modelo
#Con base en los resultados:

#Mejor Modelo: AR_1_MA_6
#AIC: 22553.76
#RMSE: 1.528931
#Log Likelihood: -11267.88
#Este modelo no solo tiene el menor AIC, sino que también presenta un RMSE bajo y el log likelihood más alto, lo que sugiere que se ajusta mejor a los datos que los otros modelos.

# Ajustar el modelo ARIMA(1, 0, 6) - AR_1_MA_6
modelo_arima_6 <- arima(na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(1, 0, 6))

# Extraer residuales
residuales <- modelo_arima_6$residuals

# Graficar los residuales
par(mfrow = c(2, 2))  # Configurar el espacio para múltiples gráficos

# 1. Gráfico de los residuales
plot(residuales, main = "Residuales del Modelo AR_1_MA_6", ylab = "Residuales", xlab = "Índice")
abline(h = 0, col = "red")

# 2. Histograma de los residuales
hist(residuales, breaks = 30, main = "Histograma de Residuales", xlab = "Residuales", border = "blue", col = "lightblue")

# 3. Gráfico de Autocorrelación (ACF)
acf(residuales, main = "Función de Autocorrelación de Residuales")

# 4. Gráfico de Autocorrelación Parcial (PACF)
pacf(residuales, main = "Función de Autocorrelación Parcial de Residuales")

# 5. Prueba de Ljung-Box para ruido blanco
ljung_box_result <- Box.test(residuales, lag = 20, type = "Ljung-Box")
print(ljung_box_result)

## 
##  Box-Ljung test
## 
## data:  residuales
## X-squared = 27.104, df = 20, p-value = 0.1324

# 6. Gráfico QQ para comprobar la normalidad de los residuales
qqnorm(residuales)
qqline(residuales, col = "red", lwd = 2)
title("Gráfico QQ de Residuales")

#Interpretación
#Hipótesis nula (H0): La hipótesis nula de la prueba de Box-Ljung es que no hay autocorrelación en los residuos (los residuos son ruido blanco).

#p-value: El p-value de 0.1324 indica que no hay suficiente evidencia para rechazar la hipótesis nula. En general, un p-value mayor que 0.05 sugiere que no hay autocorrelación significativa en los residuos.

#X-squared: El valor de 27.104 junto con los grados de libertad (20) sugiere que la estadística de prueba no es lo suficientemente grande como para proporcionar evidencia fuerte contra la hipótesis nula.

#Conclusión
#Dado que el p-value es 0.1324 (mayor que 0.05), puedes concluir que no hay evidencia significativa de autocorrelación en los residuos de tu modelo. Esto sugiere que el modelo ARIMA que has utilizado es adecuado y que los residuos se comportan como un proceso de ruido blanco, lo que es un buen signo para la validez del modelo. Esto significa que el modelo ha capturado la estructura de la serie temporal correctamente y que las predicciones realizadas por el modelo son razonablemente confiables.

# Filtrar datos para el periodo específico
datos_filtrados_periodo <- datos_filtrados[datos_filtrados$Fecha >= "2024-04-11", ]

# Estimar el modelo ARIMA(1,0,6)
modelo_arima <- arima(na.omit(datos_filtrados$Segunda_Diferencia_Precio), order = c(1, 0, 6))

# Generar pronósticos para 15 días
forecast_values <- forecast(modelo_arima, h = 15)

# Crear un data frame para las fechas de pronóstico
ultima_fecha <- max(datos_filtrados$Fecha, na.rm = TRUE)
fechas_pronostico <- seq(from = ultima_fecha + 1, 
                         by = "day", 
                         length.out = length(forecast_values$mean))

# Convertir el forecast a un data frame
forecast_df <- data.frame(Fecha = fechas_pronostico,
                           Pronóstico = forecast_values$mean)

# Graficar las segundas diferencias y el pronóstico
library(ggplot2)
library(plotly)  # Cargar el paquete plotly

# Crear el gráfico
grafico <- ggplot() +
  geom_line(data = datos_filtrados_periodo, aes(x = Fecha, y = Segunda_Diferencia_Precio, color = "Segunda Diferencia Precio")) +
  geom_line(data = forecast_df, aes(x = Fecha, y = Pronóstico, color = "Pronóstico ARIMA(1,0,6)")) +
  labs(title = "Pronóstico ARIMA(1,0,6) y Segunda Diferencia de Precio",
       x = "Año",
       y = "Segunda Diferencia Precio",
       color = "Leyenda") +
  scale_color_manual(values = c("Segunda Diferencia Precio" = "blue", "Pronóstico ARIMA(1,0,6)" = "red")) +
  theme_minimal()

# Convertir el gráfico ggplot a un gráfico interactivo de plotly
grafico_interactivo <- ggplotly(grafico)

# Mostrar el gráfico interactivo
grafico_interactivo

# TOMANDO EN CUENTA QUE EL ULTIMO PRECIO CONOCIDO ES DE 68.84 USD DEL 2024-09-19 
#AL POSICIONARNOS CON EL PUNTERO PODEMOS VER EN RELACION A ESE PRECIO CUANTO VA A VARIAR EN DOLARES
#POR EJEMPLO EL 2024-09-20 LA VARIACION ES DE -0.36637 USD DANDO UN PRECIO DE 68.47 PARA EL DIA SIGUIENTE

#Conclusiones sobre los Modelos ARIMA
#Eficiencia: ARIMA es efectivo para predecir series temporales con patrones de autocorrelación.

#Flexibilidad: Combina componentes autorregresivos y de media móvil, adaptándose a diversas situaciones.

#Estacionariedad: Requiere que la serie sea estacionaria, lo que puede necesitar diferenciación y transformaciones.

#Pronósticos Cortos: Mejor rendimiento en pronósticos a corto plazo; la precisión disminuye a largo plazo.