Master Big Data y Business Analytics - Facultad de Estudios EstadÃsticos 2022-2023
Introducción
El objetivo de este trabajo es realizar un análisis de series de tiempo sobre datos de consumo de energÃa. La base de datos proviene de la Administración de Información de Energia de los Estados Unidos de América, o E.I.A por sus siglas en inglés (https://www.eia.gov/totalenergy/data/monthly/).
En esta se pueden observar los consumos de energÃa mensuales desde el año 1943, según el tipo de consumo. Para este trabajo se ha escogido el consumo total de combustibles fósiles. La unidad de medida es en Quadrillion Btu.
El consumo de energÃa suele seguir una estacionalidad guiada por los meses de uso, sea verano o invierno, como se notará en el siguiente análisis. Ademas de ver esto en una descomposición de la serie de datos, se observarán las funciones de autocorrelaciones y autocorrelaciones parciales. Por último, se generarán modelos de predicción y se verá cual ajusta de mejor manera a la serie temporal.
Pasos previos
Eliminamos objetos previos y notación cientÃfica.
rm(list = ls())
options(scipen = 999) Cargamos librerias.
library(readxl)
library(forecast)
library(ggplot2)
library( dplyr )
library (fUnitRoots)
library(lmtest)
library(data.table)
library(stats)Lectura y visualización de los datos
Seteamos working directory y leemos data.
setwd("C:/Users/denis/OneDrive/Escritorio/Master/Modulo 8 - Time Series/Tarea")
data <- read.csv('energy.csv',header=TRUE)
data <- data[,2:3]Observamos primeros datos y resumen.
| YYYYMM | Value |
|---|---|
| 200501 | 8.015701 |
| 200502 | 7.152356 |
| 200503 | 7.551950 |
| 200504 | 6.615682 |
| 200505 | 6.708309 |
| 200506 | 6.914146 |
summary(data)## YYYYMM Value
## Min. :200501 Min. :4.972
## 1st Qu.:200906 1st Qu.:6.216
## Median :201311 Median :6.662
## Mean :201340 Mean :6.670
## 3rd Qu.:201803 3rd Qu.:6.975
## Max. :202208 Max. :8.066Contamos con información desde 2005-01 a 2022-08. No contenemos NAs y los valores se manejan dentro del rango de 4 a 8 Quad-BTUs.
Representación gráfica y descomposición de la misma
Creamos objeto time series.
#Columna 2; Comienzo 2005-01; Frequencia mensual.
energy <- ts(data[2],start=c(2005,1),frequency=12) Visualizamos serie.
autoplot(energy)+ ggtitle("Consumo de energia mensual") +
xlab("Mes") + ylab("energy") +
theme_minimal()
Descomponemos serie de tiempo con la funcion decompose(). Se pueden observar tanto los datos originales como la tendencia, la estacionalidad y el error.
Como se comentó previamente, el consumo energético tiene una componente estacional, como también puede verse en el tercer plot debajo. Asà mismo, cabe destacar la caÃda del consumo en 2020, causado principalmente por la pandemia. Debido a que la estacionalidad varÃa a lo largo del tiempo, tanto por la pandemia como mayores valores al comienzo de la serie, se escogerá el método multiplicativo.
energy_comp<- decompose(energy,type=c("multiplicative"))
autoplot(energy_comp) +
theme_minimal()
Vease como en Enero hay un 14% más consumo de energÃa que en la media anual, y en Abril un 9% menos. Luego en Junio y Julio sube arriba de la media, como tambien en Diciembre (9%).
energy_comp$seasonal[1:12]## [1] 1.1444442 1.0245120 1.0239235 0.9103513 0.9219512 0.9514534 1.0233716
## [8] 1.0316964 0.9347594 0.9530934 0.9806625 1.0997812De la misma manera, podemos ver tanto un seasonal plot del consumo de la energÃa, mostrando el consumo mensual para cada año particular, como tambien la serie desestacionalizada.
ggseasonplot(energy, year.labels=TRUE, year.labels.left=TRUE) +
ylab("Consumo") +
ggtitle("Seasonal plot: energy") + theme_minimal()
autoplot(energy, series="Datos") +
autolayer(trendcycle(energy_comp), series="Tendencia") +
autolayer(seasadj(energy_comp), series="Estacionalmente ajustada") +
xlab("Year") + ylab("Energy") +
ggtitle("Energy") +
scale_colour_manual(values=c("gray","blue","red"),
breaks=c("Datos","Estacionalmente ajustada","Tendencia")) + theme_minimal()
Separamos últimos datos para evaluar predicciones futuras.
energy_TR<-window(energy,start=c(2005,1), end=c(2020,1))
energy_Test<- window(energy,start=c(2020,2))Modelo de suavizado exponencial más adecuado
Debido a que la serie presenta estacionalidad, no tiene sentido realizar un modelo de alisado simple ni el de alisado doble de Holt
Por esto, se comenzará con el modelo Holt-Winters, con la función hw().
\(X_t = (L_t + b_t)*S_t + Z_t\)
fit_hw <- hw(energy_TR,h=12, seasonal="multiplicative")
autoplot(fit_hw) +
autolayer(fitted(fit_hw), series="Fitted") +
ylab("Energy") + xlab("Fecha") + theme_minimal()
summary(fit_hw)##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = energy_TR, h = 12, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha = 0.4305
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 7.129
## b = -0.0039
## s = 1.099 0.9786 0.953 0.9337 1.0315 1.0198
## 0.9536 0.9285 0.9167 1.0262 1.0196 1.1398
##
## sigma: 0.0275
##
## AIC AICc BIC
## 345.5990 349.3536 399.9734
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001392319 0.1803738 0.1403972 -0.02926616 2.072251 0.6467909
## ACF1
## Training set 0.03717615
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Feb 2020 6.668187 6.433577 6.902798 6.309382 7.026993
## Mar 2020 6.707867 6.450874 6.964859 6.314831 7.100903
## Apr 2020 5.988145 5.741381 6.234909 5.610752 6.365538
## May 2020 6.061636 5.795412 6.327860 5.654482 6.468791
## Jun 2020 6.222248 5.933058 6.511439 5.779970 6.664527
## Jul 2020 6.649813 6.324605 6.975021 6.152450 7.147175
## Aug 2020 6.722037 6.377721 7.066352 6.195452 7.248622
## Sep 2020 6.081316 5.756315 6.406317 5.584270 6.578362
## Oct 2020 6.203087 5.858330 6.547843 5.675827 6.730346
## Nov 2020 6.366458 5.999505 6.733410 5.805252 6.927663
## Dec 2020 7.145526 6.719430 7.571621 6.493869 7.797182
## Jan 2021 7.406042 6.950107 7.861978 6.708749 8.103336De esta manera, la expresión del modelo serÃa la siguiente:
\(L_t = 0.4305(x_t/S_{t-s}) + (1-0.4305)(L_{t-1} + b_{t-1})\)
\(b_t = 0.0001(L_t - L_{t-1}) + 0.9999b_{t-1}\)
\(S_t = 0.0001(x_t/L_t) + (1-0.0001)S_{t-s}\)
\(\hat{x}_{t+1} = (L_t + b_t)*S_{t-s+1}\)
Funciones de autocorrelación y autocorrelación parcial
Con nuestra serie de tiempo podemos observar las funciones de autocorrelacion, tanto totales como parciales, para poder encontrar patrones e importancias de los retardos del consumo de energÃa sobre el consumo actual. Con la primera podemos medir la relación lineal entre las variables de la serie separadas por k posiciones, mientras que con la segunda medimos lo mismo pero ajustando por la presencia de los demás términos de desfase más corto.
Estas pueden ser analizadas a partir de correlogramas, donde vemos la significatividad de cada retardo y los que luego serán los parámetros de nuestro modelo, como también la estacionariedad o no de nuestra serie original.
Podemos utilizar la función ggtsdisplay() para visualizar esto. Notamos que no existe estacionariedad en la serie debido a que la ACF no se reduce rápidamente, sino que muestra importancias en perÃodos estacionales como por ejemplo los 12, 24 y 36 meses de diferencia.
La PACF nos sirve para identificar el modelo autorregresivo AR(). En este caso notamos importancia hasta el periodo 12. Sin embargo, hay intervalos medios en donde no se cruza la banda de confianza.
ggtsdisplay(energy_TR,theme=theme_minimal())
Para transformar la serie a una que sea estacionaria, podemos probar generando diferencias y modelando la estacionalidad.
Primera diferencia
energy_dif<-diff(energy_TR)
ggtsdisplay(energy_dif,theme=theme_minimal())
Modelamos
modelo_dif1<-arima(energy_TR,order=c(0,1,0))
summary(modelo_dif1)##
## Call:
## arima(x = energy_TR, order = c(0, 1, 0))
##
##
## sigma^2 estimated as 0.2627: log likelihood = -135.09, aic = 272.18
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.004296411 0.5110949 0.3987326 -0.3385827 5.928287 0.9945856
## ACF1
## Training set -0.06243303Diferencia con estacionalidad anual
energy_dif12<-diff(energy_TR,12) %>% diff()
ggtsdisplay(energy_dif12,theme=theme_minimal())
PodrÃamos probar un modelo sobre esta transformación con parametros AR(2) y MA(2), siguiendo lo que muestran los gráficos de ACF y PACF.
modelo_dif12<-arima(energy_TR,order=c(2,1,2), seasonal=c(1,1,1))
summary(modelo_dif12)##
## Call:
## arima(x = energy_TR, order = c(2, 1, 2), seasonal = c(1, 1, 1))
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1
## -0.2655 -0.0831 -0.2479 -0.1871 0.2407 -0.9995
## s.e. 0.3649 0.2270 0.3610 0.3034 0.0913 0.5040
##
## sigma^2 estimated as 0.03379: log likelihood = 32.45, aic = -50.91
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.008231354 0.177116 0.1358394 0.07979956 2.014022 0.3388333
## ACF1
## Training set 0.01687959Para obtener el mejor modelo posible, sin embargo, utilizaremos la función auto_arima() del paquete forecast. Esta encuentra el mejor modelo ARIMA según los criterios de información AIC, AICc y BIC.
A este le podemos indicar la presencia de estacionalidad y luego utiliza la función nsdiffs() para determinar la D (número de diferencias estacionales a utilizar), como tambien ndiffs() para determinar la d (número de diferencias ordinarias a utilizar).
Notamos como la función escoge como mejor modelo a ARIMA(2,0,1)(2,1,1)[12]. Este no diferencia la serie para el modelo general, sino que lo hace para el modelo estacional, una vez sola. El componente AR() del modelo general tiene como parametro a p=2, y del MA() a q=1. Los mismos parametros se ven en el modelo estacional, siendo P=2 y Q=1. Este modelo muestra una mejora menor comparado con el intentado previamente.
modelo <- auto.arima(energy_TR,seasonal=TRUE)
summary(modelo)## Series: energy_TR
## ARIMA(2,0,1)(2,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 sma1
## 1.0556 -0.0877 -0.6532 0.1412 -0.2949 -0.8014
## s.e. 0.0293 0.0303 0.0615 0.0363 0.0306 0.0744
##
## sigma^2 = 0.03472: log likelihood = 39.12
## AIC=-64.24 AICc=-63.55 BIC=-42.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01782103 0.1768301 0.1351414 -0.3107178 2.005757 0.622578
## ACF1
## Training set 0.0006198699Podemos visualizar sus residuos. Vemos el comportamiento de white noise por lo que funciona correctamente. Tambien podemos confirmarlo con el test de Ljung-Box. Este tiene como hipótesis nula que los residuos estan independientemente distribuidos (son white noise). Con un P-Valor de 0.3932, no podemos rechazar H0, por lo que hay evidencia suficiente para reconocer que los errores del modelo son ruido blanco.
ggtsdisplay(residuals(modelo))
checkresiduals(modelo)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,1)(2,1,1)[12]
## Q* = 18.978, df = 18, p-value = 0.3932
##
## Model df: 6. Total lags used: 24De la misma manera, podemos observar los valores ajustados y los originales.
cbind("EnergÃa" = energy_TR, "Valores ajustados" =fitted(modelo)) %>%
autoplot() + xlab("Periodo") + ylab("") + ggtitle("Consumo de energÃa")
Expresión algebraica del modelo ajustado
Para escribir la la expresión algebráica del modelo necesitamos obtener previamente sus coeficientes.
coeftest(modelo)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.055580 0.029340 35.9778 < 2.2e-16 ***
## ar2 -0.087700 0.030278 -2.8965 0.0037733 **
## ma1 -0.653191 0.061473 -10.6256 < 2.2e-16 ***
## sar1 0.141179 0.036297 3.8896 0.0001004 ***
## sar2 -0.294897 0.030629 -9.6279 < 2.2e-16 ***
## sma1 -0.801369 0.074354 -10.7778 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1De esta manera, el modelo serÃa el siguiente:
\((1-1.055B + 0.087B^2)(1-0.141B^{12} + 0.294B^{24})(1-B^{12})X_t = (1 - 0.653B)(1 - 0.801B^{12})Z_t\)
Predicciones e intervalos de confianza
Calcularemos las predicciones para un horizonte temporal de 31 meses siguientes, hasta llegar al último mes disponible en nuestra serie de tiempo original.
autoplot(forecast(modelo),h=length(energy_Test)) + theme_minimal()
Y su summary: Notamos un AIC de -64 y un RMSE de 0.176. Asà mismo observamos las predicciones y sus intervalos de confianza en formato tabla.
summary(forecast(modelo,h=length(energy_Test)))##
## Forecast method: ARIMA(2,0,1)(2,1,1)[12]
##
## Model Information:
## Series: energy_TR
## ARIMA(2,0,1)(2,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 sma1
## 1.0556 -0.0877 -0.6532 0.1412 -0.2949 -0.8014
## s.e. 0.0293 0.0303 0.0615 0.0363 0.0306 0.0744
##
## sigma^2 = 0.03472: log likelihood = 39.12
## AIC=-64.24 AICc=-63.55 BIC=-42.34
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01782103 0.1768301 0.1351414 -0.3107178 2.005757 0.622578
## ACF1
## Training set 0.0006198699
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Feb 2020 6.671215 6.432336 6.910095 6.305881 7.036549
## Mar 2020 6.715430 6.457942 6.972919 6.321636 7.109225
## Apr 2020 5.887477 5.617699 6.157255 5.474887 6.300067
## May 2020 6.037646 5.757221 6.318072 5.608772 6.466520
## Jun 2020 6.182813 5.892857 6.472769 5.739364 6.626263
## Jul 2020 6.690408 6.391859 6.988957 6.233817 7.146999
## Aug 2020 6.710649 6.404321 7.016976 6.242162 7.179135
## Sep 2020 6.115691 5.802300 6.429082 5.636400 6.594981
## Oct 2020 6.172771 5.852948 6.492594 5.683643 6.661898
## Nov 2020 6.414770 6.089076 6.740463 5.916664 6.912875
## Dec 2020 7.199739 6.868677 7.530800 6.693424 7.706053
## Jan 2021 7.385090 7.049121 7.721059 6.871270 7.898910
## Feb 2021 6.578407 6.215744 6.941070 6.023761 7.133052
## Mar 2021 6.677975 6.305251 7.050699 6.107943 7.248008
## Apr 2021 5.957875 5.576909 6.338842 5.375238 6.540513
## May 2021 6.047415 5.659011 6.435819 5.453402 6.641428
## Jun 2021 6.253863 5.858670 6.649056 5.649468 6.858258
## Jul 2021 6.715792 6.314386 7.117198 6.101895 7.329689
## Aug 2021 6.737732 6.330631 7.144833 6.115125 7.360339
## Sep 2021 6.126007 5.713678 6.538336 5.495405 6.756610
## Oct 2021 6.222855 5.805721 6.639989 5.584904 6.860806
## Nov 2021 6.409712 5.988157 6.831267 5.765000 7.054424
## Dec 2021 7.236814 6.811189 7.662440 6.585876 7.887753
## Jan 2022 7.598524 7.169158 8.027891 6.941865 8.255184
## Feb 2022 6.625576 6.194034 7.057118 5.965590 7.285563
## Mar 2022 6.777176 6.342956 7.211395 6.113094 7.441257
## Apr 2022 6.015388 5.578630 6.452147 5.347424 6.683353
## May 2022 6.100087 5.660973 6.539200 5.428521 6.771652
## Jun 2022 6.276112 5.834819 6.717405 5.601213 6.951012
## Jul 2022 6.764918 6.321607 7.208230 6.086932 7.442905
## Aug 2022 6.795712 6.350531 7.240893 6.114866 7.476557Comparación de las predicciones obtenidas
Comparamos las predicciones contra el set test reservado previamente. Notamos que la performance empeora con el test como era de esperarse.
fit_modelo<-forecast(modelo,h=length(energy_Test),seasonal=TRUE)
accuracy(fit_modelo, x = energy_Test)## ME RMSE MAE MPE MAPE MASE
## Training set -0.01782103 0.1768301 0.1351414 -0.3107178 2.005757 0.622578
## Test set -0.18458703 0.3291021 0.2174863 -3.2363344 3.721505 1.001930
## ACF1 Theil's U
## Training set 0.0006198699 NA
## Test set 0.6819457206 0.7198352Sin embargo, las predicciones son acertadas, sin contar el año 2020 donde la pandemia tiene un peso negativo.
En rojo, los valores reales, en negro las predicciones.
fit_modelo %>%
autoplot() +
geom_line(
aes(
x = as.numeric(time(energy_Test)),
y = as.numeric(energy_Test)
),
col = "red"
)
Podemos aplicar el modelo al set de testing y ver que los fitted values son cercanos a los originales.
modelo_test <- Arima(energy_Test, model = modelo)
one_step_forecasts <- fitted(modelo_test)
ggplot() +
geom_line(
aes(
x = as.numeric(time(one_step_forecasts)),
y = as.numeric(one_step_forecasts)
),
col = "black"
) +
geom_line(
aes(
x = as.numeric(time(energy_Test)),
y = as.numeric(energy_Test)
),
col = "red"
)