Análisis del consumo de gas anual en México

El GLP es uno de los energéticos más utilizados en México, al punto de ser identificada como la nación con el mayor consumo anual. El uso del gas LP es fundamental para los consumidores mexicanos. A partir de su uso en los hogares se puede cocinar alimentos y calentar agua.

Análisis de series de tiempo (TSA)

A continuacion se mostrara un análisis de series de tiempo el cual nos permite analizar el comportamiento del consumo anual del GLP en México de 1966 hasta 2019.

Datos y Paquetes

library(pacman)
p_load("tidyverse", "lubridate", "forecast", "TTR", "MLmetrics", "tseries", "fpp", "TSstudio", "xfun", "plotly")

dataGasMexico <- read_csv("annual-change-gas_mexico.csv")

## Rows: 54 Columns: 4

## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (2): Entity, Code
## dbl (2): Year, Gas

## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

Importando y graficando datos

La siguiente grafica muestra la serie de tiempo del consumo anual que tiene el Gas LP en México, empezando del año 1966 hasta el 2019

datos<-ts(dataGasMexico$Gas, frequency = 2 ,start = c(1966, 1),end = c(2019,1))
e <- plot_ly(data = dataGasMexico, x = dataGasMexico$Year , y = dataGasMexico$Gas, mode = 'lines', type="scatter", name="Gas") %>%  layout(title = 'Consumo de gas en México' ,
         yaxis = list(zeroline = FALSE, title = 'TWh'),
         xaxis = list(zeroline = FALSE, title = 'Año'), legend = list(font = list(size = 15)))
e

Buscando datos Faltantes

anyNA(datos)

## [1] FALSE

Esto significa no hay valores faltante

Descomposición de la serie de tiempo

La descomposición de series de tiempo deconstruye una serie de tiempo en varios componentes, cada uno de los cuales representa una de las categorías subyacentes de patrones

datos_dc <- datos %>%
  decompose(type="multiplicative") %>%
  autoplot()
datos_dc

1. Usando un modelo de suavizado exponencial

• Validación cruzada

La validación cruzada o cross-validation es una técnica utilizada para evaluar los resultados de un análisis estadístico y garantizar que son independientes de la partición entre datos de entrenamiento y prueba.

test_datos <- tail(datos, 11)  #20% para pruebas
train_datos <- head(datos, length(datos)-11) #80% para entrenamiento

• Holt-Winters

datos_ses <- HoltWinters(train_datos, seasonal = "multiplicative", )

## Warning in HoltWinters(train_datos, seasonal = "multiplicative", ): optimization
## difficulties: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH

datos_ses

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = train_datos, seasonal = "multiplicative")
## 
## Smoothing parameters:
##  alpha: 0.3
##  beta : 0.1
##  gamma: 0.1
## 
## Coefficients:
##          [,1]
## a  22.3379643
## b   0.6786307
## s1  1.7578945
## s2  0.7117851

Forecasting

datos_forecast <- forecast(object = datos_ses, h = 11)

Visualizando resultados

datos %>%
  autoplot(series = "Datos Reales") +
  autolayer(datos_forecast$fitted, series = "Datos Entrenamiento") +
  autolayer(datos_forecast$mean , series = "Datos de Prueba") +
  theme_minimal()

## Warning: Removed 2 row(s) containing missing values (geom_path).

Evaluar la precisión del modelo

eval_ses <- accuracy(datos_forecast, test_datos)
eval_ses

##                     ME     RMSE      MAE        MPE     MAPE     MASE
## Training set  1.232047 21.83827 16.17108 -508.42008 657.0695 1.047354
## Test set     -6.110332 26.56957 23.67103  -89.67134 128.8620 1.533103
##                    ACF1 Theil's U
## Training set  0.4178301        NA
## Test set     -0.6841288 0.9019048

Viendo la prueba con los datos, el modelo tiene un promedio de error del 128%, por lo que este no es capaz de predecir el consumo del gas

Utilizando el 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 (P,D,Q)

Empezaremos haciendo una prueba al conjunto de datos

adf.test(datos)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  datos
## Dickey-Fuller = -3.409, Lag order = 4, p-value = 0.0567
## alternative hypothesis: stationary

De acuerdo con el valor de P, los datos no son estacionarios

datos_diff <- datos

datos_diff %>%
  autoplot()

datos_auto <- auto.arima(y = train_datos, seasonal = T)
summary(datos_auto)

## Series: train_datos 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     mean
##       0.7965  -0.5511  13.5557
## s.e.  0.1243   0.1645   3.5062
## 
## sigma^2 = 263.9:  log likelihood = -402.43
## AIC=812.87   AICc=813.31   BIC=823.13
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1897213 15.99028 12.08065 -840.2523 937.6974 0.7824285
##                     ACF1
## Training set -0.03206988

train_datos %>%
  autoplot(series= "Datos Reales") +
  autolayer (datos_auto$fitted, series = "SARIMA auto") +
  theme_minimal()

tsdisplay(datos_diff)

xfun::embed_file("Caso de Estudio 5.Rmd")