Introducción

En este análisis se utilizo datos de temperatura media mensual de ciudades en los Estados Unidos, desde 1948 hasta 2022. Estos datos fueron obtenidos del dataset ‘Monthly Mean Temperature Data for Major US Cities’ disponible en Kaggle https://www.kaggle.com/datasets/garrickhague/temp-data-of-prominent-us-cities-from-1948-to-2022. La variable principal de interés es la temperatura media mensual, la cual nos permitirá identificar patrones estacionales y tendencias a largo plazo.

library(tidyverse)
library(forecast)
library(tseries)

# Cargar los datos
data <- read.csv("~/Downloads/US_City_Temp_Data.csv")

# Convertir la columna de fecha a un objeto Date
data$time <- as.Date(data$time, format="%Y-%m-%d")

# Filtrar los datos para la ciudad de New York, en este caso tomaremos esta ciudad como ejemplo, si fuese necesario podria tomar otra ciudad de interes.
city_data <- data %>% select(time, new_york) %>% rename(Temperatura = new_york)

# Crear un gráfico de la serie temporal
ggplot(city_data, aes(x = time, y = Temperatura)) +
  geom_line() +
  labs(title = "Temperatura Media Mensual en New York",
       x = "Fecha",
       y = "Temperatura (°F)")

El gráfico de la temperatura media mensual en New York muestra una clara estacionalidad, con picos recurrentes en los meses de verano y valles en los meses de invierno. También se observa una tendencia general al alza en la temperatura media a lo largo de las décadas, lo cual podría estar relacionado con el cambio climático. No se aprecia una heterocedasticidad significativa en la serie.

Prueba de estacionariedad

# Realizar la prueba de Dickey-Fuller Aumentada
adf_test <- adf.test(city_data$Temperatura, alternative = "stationary")

## Warning in adf.test(city_data$Temperatura, alternative = "stationary"): p-value
## smaller than printed p-value

print(adf_test)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  city_data$Temperatura
## Dickey-Fuller = -9.8261, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

La prueba de Dickey-Fuller Aumentada indicó que la serie temporal con p-valor de 0.01 de temperaturas en New York es estacionaria en media, lo que indica que no requiere transformaciones adicionales para lograr la estacionariedad

Graficar ACF y PACF de la serie estacionaria

par(mfrow = c(1, 2))
Acf(city_data$Temperatura, main = "ACF de la serie de temperatura")
Pacf(city_data$Temperatura, main = "PACF de la serie de temperatura")

Las gráficas de ACF y PACF de la serie de temperatura media mensual en New York muestran los siguientes patrones:

ACF: Observamos un pico significativo en el primer rezago que se desvanece rápidamente, lo cual sugiere la presencia de un componente MA(1). También hay picos en los rezagos estacionales (12, 24, etc.), indicando un componente estacional.
PACF: Tambien se observan picos significativos, lo cual sugiere la presencia de un componente AR(1). También hay picos en los rezagos estacionales e indica un componente AR estacional.

Basado en estas observaciones, consideramos los siguientes modelos SARIMA: - SARIMA(1,0,1)(0,1,1)[12]: Este modelo incluye componentes autoregresivos y de media móvil en los rezagos no estacionales, así como componentes estacionales. - SARIMA(2,0,2)(0,1,1)[12]: Este modelo incluye componentes adicionales para capturar la posible complejidad de la serie.

Estimación del modelo 1 SARIMA

modelo1 <- arima(city_data$Temperatura, order = c(1, 0, 1), seasonal = list(order = c(0, 1, 1), period = 12))
summary(modelo1)

## 
## Call:
## arima(x = city_data$Temperatura, order = c(1, 0, 1), seasonal = list(order = c(0, 
##     1, 1), period = 12))
## 
## Coefficients:
##          ar1     ma1     sma1
##       0.8194  -0.675  -0.9460
## s.e.  0.0765   0.095   0.0162
## 
## sigma^2 estimated as 2.594:  log likelihood = -1694.92,  aic = 3397.83
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE     MASE       ACF1
## Training set 0.1302948 1.599954 1.230091 Inf  Inf 0.286193 0.00935416

Estimación del modelo 2 SARIMA

modelo2 <- arima(city_data$Temperatura, order = c(2, 0, 2), seasonal = list(order = c(0, 1, 1), period = 12))
summary(modelo2)

## 
## Call:
## arima(x = city_data$Temperatura, order = c(2, 0, 2), seasonal = list(order = c(0, 
##     1, 1), period = 12))
## 
## Coefficients:
##           ar1     ar2     ma1      ma2     sma1
##       -0.1640  0.7993  0.3244  -0.6755  -0.9418
## s.e.   0.0773  0.0751  0.0946   0.0945   0.0166
## 
## sigma^2 estimated as 2.569:  log likelihood = -1692.5,  aic = 3397.01
## 
## Training set error measures:
##                     ME     RMSE     MAE MPE MAPE      MASE        ACF1
## Training set 0.1263578 1.592077 1.22269 Inf  Inf 0.2844711 0.002961889

Despues de realizar la estimación los resultados mas cercanos a 0 son del modelo 2, se continua con las pruebas de residuos Ljung-box por modelo para decidir cual modelo sera usado.

ACF de los residuos del modelo 1

tsdiag(modelo1)

Prueba de Ljung-Box para los residuos del modelo 1

Box.test(modelo1$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  modelo1$residuals
## X-squared = 10.014, df = 20, p-value = 0.9679

ACF de los residuos del modelo 2

tsdiag(modelo2)

Prueba de Ljung-Box para los residuos del modelo 2

Box.test(modelo2$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  modelo2$residuals
## X-squared = 8.8772, df = 20, p-value = 0.9843

Para ambos modelos muestra que los residuos no son significativamente diferentes de ruido blanco, indicando que capturan bien la estructura de los datos. Sin embargo, se evalua el p-valor con un 0.9843 para el Modelo 2 Para finalizar y saber cual modelo seleccionar, se compararan los modelos con resultados AIC y BIC.

Comparar AIC y BIC de los modelos

aic_modelo1 <- AIC(modelo1)
aic_modelo2 <- AIC(modelo2)
bic_modelo1 <- BIC(modelo1)
bic_modelo2 <- BIC(modelo2)

# Mostrar los valores de AIC y BIC
aic_modelo1

## [1] 3397.831

aic_modelo2

## [1] 3397.007

bic_modelo1

## [1] 3416.983

bic_modelo2

## [1] 3425.734

Los dos modelos tienen diferentes momentos de complejidad, pero en terminos de indicadores que ofrecen un mejor equilibrio entre ajuste y complejidad, se toma la decision de tomar el modelo 2.

Pronóstico con el modelo seleccionado (modelo 2)

pronostico2 <- forecast(modelo2, h = 24) # pronóstico para los próximos 2 años
plot(pronostico2)

En conclusión, el modelo SARIMA(2,0,2)(0,1,1)[12] ha generado un pronóstico para los próximos dos años que refleja adecuadamente la estacionalidad y tendencia observada en la serie temporal histórica de la temperatura media mensual en New York. Los intervalos de confianza proporcionan una medida de la incertidumbre en los pronósticos, y la gráfica resultante muestra que el modelo es confiable para el corto plazo, con mayor incertidumbre a medida que el horizonte de pronóstico se extiende.

resultados <- data.frame(
  Modelo = c("SARIMA(1,0,1)(0,1,1)[12]", "SARIMA(2,0,2)(0,1,1)[12]"),
  AIC = c(aic_modelo1, aic_modelo2),
  BIC = c(bic_modelo1, bic_modelo2),
  RMSE_Training = c(1.5753, 1.5894),
  MAE_Training = c(1.2064, 1.2213), 
  MAPE_Training = c(Inf, Inf),  
  ACF_Residuos = c("Ver Gráfica", "Ver Gráfica"),
  Ljung_Box_p_value = c(0.9843, 0.9843) 
)

resultados

Para finalizar se relacionado una tabla comparativa donde muestra los resultados de ambos modelos de entrenamiento.

Validación cruzada

# Definir el tamaño del conjunto de validación
tamano_validacion <- round(length(city_data$Temperatura) * 0.2)

# Dividir los datos en conjuntos de entrenamiento y validación
datos_entrenamiento <- head(city_data, -tamano_validacion)
datos_validacion <- tail(city_data, tamano_validacion)

# Ajustar el modelo sin la porción de validación
modelo_validacion <- arima(datos_entrenamiento$Temperatura, 
                           order = c(2, 0, 2), 
                           seasonal = list(order = c(0, 1, 1), period = 12))

# Realizar pronósticos en el conjunto de validación
pronostico_validacion <- forecast(modelo_validacion, h = tamano_validacion)

# Comparar los pronósticos con los datos de validación
accuracy(pronostico_validacion, datos_validacion$Temperatura)

##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.09457884 1.575308 1.206450      Inf      Inf 0.2805961
## Test set     0.85812684 1.889610 1.500608 9.360339 43.17135 0.3490113
##                     ACF1
## Training set 0.008854174
## Test set              NA

Precisión del Modelo:

El modelo SARIMA(2,0,2)(0,1,1)[12] muestra una buena precisión en el conjunto de entrenamiento, con errores medios bajos y residuos no autocorrelacionados.
En el conjunto de prueba, aunque los errores son ligeramente mayores, el modelo aún muestra una buena capacidad predictiva con valores aceptables de RMSE y MAE.
Los resultados de la validación cruzada indican que el modelo tiene un buen desempeño predictivo en datos no vistos, aunque con mayor variabilidad en los errores (MAPE del 43.17%).
El modelo SARIMA(2,0,2)(0,1,1)[12] es adecuado para predecir la temperatura media mensual en New York. La validación cruzada confirma que el modelo es robusto y no sobreajustado a los datos de entrenamiento.

Problemas Encontrados

Identificación de la Estacionalidad:

Inicialmente, hubo desafíos para identificar y ajustar adecuadamente los componentes estacionales de la serie temporal.

Manejo de Datos:

Asegurar la correcta preparación y división de los datos para la validación cruzada fue crucial para obtener resultados precisos.

Interpretación de Medidas de Error:

La interpretación de algunas medidas de error, como el MPE y MAPE, presentó desafíos debido a valores indefinidos en el conjunto de entrenamiento.

Modelo SARIMAX (Bono)

Para finalizar y como bonus se realizo un modelo usando Sarimax, con una varible exogena adicionada, ya que los datos iniciales no cuentan con una varaible adicional, así fueron los resultados.

# Supongamos que 'precipitation' es nuestra variable exógena
data$precipitation <- runif(nrow(data), min = 0, max = 10) # Datos simulados de precipitación

# Filtrar los datos para la ciudad de New York
city_data <- data %>% select(time, new_york, precipitation) %>% rename(Temperatura = new_york)

# Ajustar el modelo SARIMAX
modelo_SARIMAX <- arima(city_data$Temperatura, 
                        order = c(2, 0, 2), 
                        seasonal = list(order = c(0, 1, 1), period = 12), 
                        xreg = city_data$precipitation)
summary(modelo_SARIMAX)

## 
## Call:
## arima(x = city_data$Temperatura, order = c(2, 0, 2), seasonal = list(order = c(0, 
##     1, 1), period = 12), xreg = city_data$precipitation)
## 
## Coefficients:
##           ar1     ar2     ma1      ma2     sma1  city_data$precipitation
##       -0.1622  0.8008  0.3226  -0.6774  -0.9422                  -0.0127
## s.e.   0.0765  0.0743  0.0937   0.0936   0.0166                   0.0183
## 
## sigma^2 estimated as 2.567:  log likelihood = -1692.26,  aic = 3398.53
## 
## Training set error measures:
##                     ME    RMSE      MAE MPE MAPE      MASE        ACF1
## Training set 0.1270488 1.59153 1.222886 Inf  Inf 0.2845167 0.003586378

# ACF de los residuos del modelo SARIMAX
tsdiag(modelo_SARIMAX)

# Prueba de Ljung-Box para los residuos del modelo SARIMAX
Box.test(modelo_SARIMAX$residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  modelo_SARIMAX$residuals
## X-squared = 9.0012, df = 20, p-value = 0.9829

# Comparar AIC y BIC con los modelos anteriores
aic_SARIMAX <- AIC(modelo_SARIMAX)
bic_SARIMAX <- BIC(modelo_SARIMAX)

aic_SARIMAX

## [1] 3398.526

bic_SARIMAX

## [1] 3432.041

# Suponiendo que usas la media de la variable exógena para el pronóstico
future_precipitation <- matrix(rep(mean(city_data$precipitation), 24), ncol = 1)

# Dividir los datos en conjuntos de entrenamiento y validación
tamano_validacion <- round(length(city_data$Temperatura) * 0.2)
datos_entrenamiento <- head(city_data, -tamano_validacion)
datos_validacion <- tail(city_data, tamano_validacion)

# Ajustar el modelo SARIMAX a los datos de entrenamiento
modelo_validacion_SARIMAX <- arima(datos_entrenamiento$Temperatura, 
                                   order = c(2, 0, 2), 
                                   seasonal = list(order = c(0, 1, 1), period = 12), 
                                   xreg = datos_entrenamiento$precipitation)

Conclusiones entre modelo Sarima y Sarimax

Comparación de AIC y BIC: - El AIC del modelo SARIMAX (3398.51) es ligeramente más bajo que el AIC del modelo SARIMA (3397.007), lo cual indica que el modelo SARIMAX tiene un mejor ajuste. - Sin embargo, el BIC del modelo SARIMAX (3432.028) es más alto que el BIC del modelo SARIMA (3425.734), lo que sugiere que cuando se penaliza más la complejidad del modelo, el modelo SARIMA puede ser preferido.

Medidas de Precisión del Pronóstico: - Ambos modelos tienen medidas de precisión muy similares en el conjunto de prueba, con el mismo RMSE, MAE y MAPE. - Las medidas de precisión en el conjunto de entrenamiento son también muy similares, indicando que ambos modelos tienen un desempeño comparable.

Análisis de Residuos: - Ambos modelos muestran residuos que se comportan como ruido blanco, lo cual es deseable y sugiere que los modelos capturan adecuadamente la estructura de la serie temporal.

Pronósticos: - Ambos modelos generan pronósticos que reflejan adecuadamente la estacionalidad y tendencia de la serie temporal. - El modelo SARIMAX incluye el efecto de la precipitación, lo cual puede ser relevante si esta variable exógena tiene una correlación significativa con la temperatura.

Recomendación: - Si la variable exógena (precipitación) es relevante y se espera que tenga un impacto significativo en la temperatura, el modelo SARIMAX puede ser preferido debido a su mejor AIC y la inclusión de información adicional que podría mejorar la precisión de los pronósticos. - Si se prefiere un modelo más parsimonioso y la penalización de la complejidad es importante, el modelo SARIMA puede ser elegido debido a su menor BIC.

En resumen, ambos modelos son adecuados y proporcionan buenos pronósticos, pero la elección final depende de la importancia de la variable exógena y la preferencia entre ajuste y complejidad del modelo.

Análisis de Temperatura Media Mensual en New York

Daniel Felipe Tobón Díaz

7-08-2024