Caso de estudio: Temperatura del aire en Obregón
El estado de Sonora tiene una clima prominentemente seco, con el 46.5% del territorio presenta dentro de estas condiciones atmosféreicas. Puesto que se requiere analizar datos cíclicos, de determinó que las condiciones climatológicas podrían servir para el presente caso de estudio.
library(pacman)
p_load("tidyverse", "lubridate", "forecast", "TTR", "MLmetrics", "tseries", "fpp", "TSstudio", "prettydoc", "DT", "xfun", "readr", "ggplot2", "scales", "base64enc", "mime", "dplyr", "ggplot2")Importando datos
Se utilizó una base de datos de distintos parámetros de condiciones atmosféricas en periodos horarios por un periodo de 4 días en Ciudad Obregón.
Se obtuvieron los datos de CONAGUA: https://smn.conagua.gob.mx/es/observando-el-tiempo/estaciones-meteorologicas-automaticas-ema-s
tempOBR <- read.csv("Estacion_CD_OBREGON_1_semana.csv")
datatable(tempOBR)Pasar a linea de tiempo
TSTemp <- ts(data=tempOBR$Temperatura.del.Aire...C., frequency=12, start=c(2000))La función ts no permitió utilizar periodos horarios, por lo tanto, para poder continuar con el análisis cíclico se tomaron los datos como promedios anuales.
autoplot(TSTemp)
Descomponer la serie de tiempo
Temp_dc <- TSTemp %>%
decompose(type = "multiplicative") %>%
autoplot()
Temp_dc
División de datos para validación cruzada
test_temp <- tail(TSTemp, 50) #20% para pruebas
train_temp <- head(TSTemp, length(sales)-50) #80% para entrenamientoHolt-Winters
temp_ses <- HoltWinters(train_temp, seasonal = "multiplicative",)
temp_ses## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = train_temp, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.7933074
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 14.5103317
## b -0.4509907
## s1 0.7384315
## s2 0.7023152
## s3 0.7041197
## s4 0.7778658
## s5 0.9079988
## s6 1.0657861
## s7 1.1998142
## s8 1.2534435
## s9 1.2475408
## s10 1.1332963
## s11 1.0963585
## s12 1.0199629Forecasting
tempo_forecast <- forecast(object = temp_ses, h=40)Visualizando datos
TSTemp %>%
autoplot(series = "actual") +
autolayer(tempo_forecast$fitted, series = "train") +
autolayer(tempo_forecast$mean, series = "test") +
theme_minimal()## Warning: Removed 12 row(s) containing missing values (geom_path).
Se oberva que los valores de prueba van bajando mucho más que los valores actuales que se pueden observar en el conjunto de datos.
Evaluar la precisión del modelo
eval_temp <- accuracy(tempo_forecast, test_temp)
eval_temp## ME RMSE MAE MPE MAPE MASE
## Training set -0.1018421 4.011674 3.032332 0.1556122 17.55449 0.2258226
## Test set 12.2497894 15.146499 12.372306 64.2256435 65.28681 0.9213852
## ACF1 Theil's U
## Training set 0.4859783 NA
## Test set 0.8794355 6.140938Se tiene un error medio de 36.84%, lo cual es alto si se busca un modelo preciso
Modelo SARIMA o ARIMA
Los modelos SARIMA captan el comportamiento puramente estacional de una serie, en forma similar, como hemos visto, se realiza para la componente regular o no estacional. Una serie con influencia solamente por la componente estacional puede ser descrito por un modelo SARIMA
Empezamos haciendo una prueba al conjunto de datos:
Esta prueba (dickey-fuller) sirve para en funcion del valor de P, determinar si los datos SON o NO SON estacionarios
adf.test(TSTemp)## Warning in adf.test(TSTemp): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: TSTemp
## Dickey-Fuller = -5.2512, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryEl p-value es de 0.01, por lo tanto los datos son NO estacionarios.
Análisis diferencial entre ciclos (estaciones)
temp_diff <- TSTemp %>%
diff(lag=12) %>%
diff(lag=1)
temp_diff %>%
autoplot()
Utilizando un modelo SARIMA
temp_auto <- auto.arima(y= train_temp, seasonal = T)
summary(temp_auto)## Series: train_temp
## ARIMA(3,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 sar1 mean
## 2.2059 -1.5472 0.2919 -0.9094 0.3227 17.2708
## s.e. 0.1277 0.2467 0.1289 0.0541 0.1275 0.4695
##
## sigma^2 = 2.277: log likelihood = -146.49
## AIC=306.98 AICc=308.53 BIC=323.65
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04386786 1.451393 1.097834 0.1683268 8.746687 0.08175741
## ACF1
## Training set 0.00121366Gráfica de datos reales vs. Arima
train_temp %>%
autoplot(series="actual") +
autolayer(temp_auto$fitted, series ="SARIMA auto") +
theme_minimal()
Conclusión
El modelo SARIMA creo un modelo más apegado a los valores actuales, sin embargo, se espera cierta variabilidad del modelo (este tiene un error medio de 36.84%) por la naturaleza de los cambios cíclicos meteorológicos. A la fecha de realización de este trabajo, Ciudad Obregón se encuentra en una transición de estaciones, lo cual afecta la estabilidad de los datos, aunque sea solo por cuatro días.
Descarga este código
xfun::embed_file("Caso de Estudio 5 Desastres Naturales.Rmd")Descarga los datos utilizados
xfun::embed_file("Estacion_CD_OBREGON_1_semana.csv")