Librerías e Importación de datos

Utilizaremos diferentes librerías disponibles en rstudio, las cuales nos ayudarán en la manipulación y transformación de datos así como los ajustes de series temporales (SARIMA) y su validación de los supuestos estadísticos necesarios.

library(pacman)
## Warning: package 'pacman' was built under R version 4.5.3
p_load(readr, dplyr, tidyr, lubridate, ggplot2, forecast, tseries, performance)
datos = read_delim("Macusani-2017-2025.csv", 
                    delim = ";", 
                    locale = locale(decimal_mark = "."),
                    show_col_types = FALSE)
## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
##   dat <- vroom(...)
##   problems(dat)

Limpieza y transformación de datos

Se transformaron los datos de acuerdo al formato de fechas requerido para ser manejable luego para las series de tiempo, además se seleccionaron las variables necesarias. Por consiguiente, se halló un error o datos faltantes en la toma de los mismos en el año 2020 (Desde el 17 de Marzo hasta el 10 de Junio del mismo año)presumiblemente debido a la pandemia, de esta manera se evitarán sesgos en las estimaciones del modelo.

datos = datos %>%
  mutate(Fecha = make_date(year = Año, month = Mes, day = Día)) %>%
  select(Fecha, Temp_Max = `Temp. Máx (°C)`, 
         Temp_Min = `Temp. Mín (°C)`, 
         Humedad = `Humedad (%)`, 
         Precipitacion = `Precipitación (mm)`)
datos = datos %>%
  distinct(Fecha, .keep_all = TRUE)
datos_c = datos %>%
  filter(!(Fecha >= "2020-03-17" & Fecha <= "2020-06-10"))

Análisis Exploratorio

ggplot(datos_c, aes(x = Fecha, y = Temp_Min)) +
  geom_line(color = "steelblue", alpha = 0.7) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Comportamiento histórico de la Temperatura Mínima",
       subtitle = "Macusani (2017-2025). Línea roja = umbral de helada (0°C)",
       x = "Fecha", y = "Temp. Mínima (°C)") +
  theme_minimal()

df_heladas = datos_c %>%
  mutate(Helada = ifelse(Temp_Min <= 0, 1, 0),
         Mes = month(Fecha, label = TRUE),
         Año = year(Fecha))


heladas_mensual = df_heladas %>%
  group_by(Año, Mes) %>%
  summarise(Dias_Helada = sum(Helada, na.rm = TRUE), .groups = 'drop') %>%
  group_by(Mes) %>%
  summarise(Promedio_Heladas = mean(Dias_Helada, na.rm = TRUE))


ggplot(heladas_mensual, aes(x = Mes, y = Promedio_Heladas)) +
  geom_col(fill = "steelblue") +
  labs(title = "Promedio histórico de días con helada por mes",
       x = "Mes", y = "Días promedio con helada (Temp. Min <= 0°C)") +
  theme_minimal()

## Modificación temporal Se modifica y agrega el parámetro como promedio de las temperaturas mínimas comparándolas con las mínima en cada mes. Luego se modela la serie de tiempo (ts).

df_mensual = datos_c %>%
  mutate(AnioMes = floor_date(Fecha, "month")) %>%
  group_by(AnioMes) %>%
  summarise(Temp_Min_Promedio = mean(Temp_Min, na.rm = TRUE),
            Temp_Min_Absoluta = min(Temp_Min, na.rm = TRUE),
            .groups = 'drop') %>%
  filter(AnioMes >= "2017-01-01")

serie_mensual = ts(df_mensual$Temp_Min_Promedio, 
                    start = c(2017, 1), 
                    frequency = 12)

autoplot(serie_mensual) +
  geom_line(color = "darkblue") +
  labs(title = "Temperatura Mínima Promedio Mensual",
       subtitle = "Serie con estacionalidad anual (2017-2025)",
       x = "Año", y = "Temp. Mín. Promedio (°C)") +
  theme_minimal()

Luego, se observa la estacionariedad en 12 meses (dado que la estimación se dará anual) reduciendo así el ruido estadístico para las estimaciones posteriores.

descomposicion = decompose(serie_mensual, type = "multiplicative")
autoplot(descomposicion)

Rigor estadístico del modelo(pruebas de estacionariedad con Dickey Fuller)

estacio_test = adf.test(serie_mensual, alternative = "stationary")
## Warning in adf.test(serie_mensual, alternative = "stationary"): p-value smaller
## than printed p-value
estacio_test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie_mensual
## Dickey-Fuller = -7.2068, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Esta prueba confirma la estacionariedad de la serie; sin embargo, para modelar adecuadamente la estructura de tiempo anual se utilizará auto.arima.

modelo_sarima = auto.arima(serie_mensual, 
                            seasonal = TRUE, 
                            stepwise = FALSE,
                            approximation = FALSE)

summary(modelo_sarima)
## Series: serie_mensual 
## ARIMA(0,0,2)(1,1,0)[12] 
## 
## Coefficients:
##          ma1     ma2     sar1
##       0.5636  0.4310  -0.5528
## s.e.  0.0962  0.0885   0.0918
## 
## sigma^2 = 3.066:  log likelihood = -186.99
## AIC=381.98   AICc=382.43   BIC=392.16
## 
## Training set error measures:
##                       ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.03122757 1.622318 1.089705 64.99387 165.8544 0.6318011
##                    ACF1
## Training set 0.01690472
checkresiduals(modelo_sarima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,2)(1,1,0)[12]
## Q* = 11.713, df = 18, p-value = 0.8617
## 
## Model df: 3.   Total lags used: 21

Obtenemos una validez del modelo (regular: (0, 0, 2), estacional: (1, 1, 0)). Esto significa que la serie es correcta y que error de los dos meses previos influye en el mes actual, sin embargo, el modelo corrige el patrón que se obtiene que Julio siempre tenga las temperaturas más frías para que no sean iguales todos los años, ya que estos varían ligeramente.

Pronóstico

forecast_2027 = forecast(modelo_sarima, h = 12)
df_pronostico = data.frame(
  Mes = month.name,
  Pronostico = round(as.numeric(forecast_2027$mean), 2),
  LI_80 = round(as.numeric(forecast_2027$lower[, 1]), 2),
  LS_80 = round(as.numeric(forecast_2027$upper[, 1]), 2)
)

print(df_pronostico)
##          Mes Pronostico  LI_80 LS_80
## 1    January       1.37  -0.87  3.61
## 2   February       0.94  -1.64  3.51
## 3      March       0.14  -2.61  2.89
## 4      April      -1.28  -4.03  1.47
## 5        May      -3.88  -6.63 -1.12
## 6       June      -7.37 -10.12 -4.62
## 7       July      -9.81 -12.56 -7.06
## 8     August      -7.18  -9.93 -4.43
## 9  September      -4.33  -7.08 -1.58
## 10   October      -1.42  -4.17  1.34
## 11  November       0.47  -2.28  3.22
## 12  December       0.80  -1.95  3.55
autoplot(forecast_2027) +
  labs(title = "Pronóstico de Temperatura Mínima Mensual para 2027",
       subtitle = "Modelo ARIMA(0,0,2)(1,1,0)[12] con IC al 80%",
       x = "Año", y = "Temp. Mín. Promedio (°C)") +
  theme_minimal() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red")