Los precios de la gasolina siempre son un problema en Turquía, empecé a investigar los precios de los últimos siete años y cómo evolucionaron. He usado gasolina sin plomo 95 de Petrol Ofisi, que es una empresa de distribución de combustible en Turquía. Organicé los precios mensualmente entre 2013 y 2020 como la moneda turca, Lira. El conjunto de datos los podemos encontrar aquí 

library(readxl)
gasoline_df <- read_excel("gasoline_df.xlsx")
## New names:
## * `` -> ...1
head(gasoline_df)
## # A tibble: 6 x 3
##   ...1  date                gasoline
##   <chr> <dttm>                 <dbl>
## 1 1     2013-01-01 00:00:00     4.67
## 2 2     2013-02-01 00:00:00     4.85
## 3 3     2013-03-01 00:00:00     4.75
## 4 4     2013-04-01 00:00:00     4.61
## 5 5     2013-05-01 00:00:00     4.64
## 6 6     2013-06-01 00:00:00     4.72

En primer lugar, quería ver cómo evolucionaban los precios de la gasolina durante el período y analizar las líneas de tendencia ajustadas.

Debido a que la variable de respuesta se mide por una logarítmica en el modelo exponencial, tenemos que obtener el exponente de la misma y crear un marco de datos diferente para la línea de tendencia exponencial, por lo tanto:

#exponential trend model
exp.model <- lm(log(gasoline)~date,data = gasoline_df)
exp.model.df <- data.frame(x=gasoline_df$date,
                           y=exp(fitted(exp.model)))

Ahora, creamos las lineas de tendencia con ggplot usando la librería tidyverse:

#Trend lines 
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
ggplot(gasoline_df, aes(x = date, y = gasoline)) + geom_point() +
  stat_smooth(method = 'lm', aes(colour = 'linear'), se = FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,2), aes(colour = 'quadratic'), se= FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,3), aes(colour = 'cubic'), se = FALSE)+
  stat_smooth(data=exp.model.df, method = 'loess',aes(x,y,colour = 'exponential'), se = FALSE)

Como vemos arriba, los modelos cúbicos y exponenciales casi se superponen y parecen ajustarse mejor a los datos.

Los modelos de tendencia de pronóstico

La tendencia lineal ; y_{t} , el valor de la serie en un momento dado, llamado t.
Lo podemos crear de este modo:

model_linear <- lm (data = gasoline_df, gasoline ~ date)
model_linear
## 
## Call:
## lm(formula = gasoline ~ date, data = gasoline_df)
## 
## Coefficients:
## (Intercept)         date  
##  -9.941e+00    1.037e-08
Para comparar los modelos, tenemos que extraer los coeficientes de determinación ajustados, que se utilizan para comparar modelos de regresión con un número diferente de variables explicativas de cada modelo de tendencia.

Equilibrado

n : tamaño de la muestra k : número de variables R^2 : coeficiente de determinación que es el cuadrado de correlación de los valores de observación con los valores predichos. 

adj_r_squared_linear <- summary(model_linear) %>% 
  .$adj.r.squared %>% 
  round(4)

adj_r_squared_linear
## [1] 0.6064
La tendencia exponencial

Esta permite que la serie aumente a un ritmo creciente en cada período, es un logaritmo natural de la variable de respuesta.

model_exponential <- lm(data = gasoline_df,log(gasoline)~date)
summary(model_exponential)
## 
## Call:
## lm(formula = log(gasoline) ~ date, data = gasoline_df)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.216180 -0.083077  0.009544  0.087953  0.151469 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.076e+00  2.496e-01  -4.311  4.5e-05 ***
## date         1.859e-09  1.701e-10  10.928  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09939 on 82 degrees of freedom
## Multiple R-squared:  0.5929, Adjusted R-squared:  0.5879 
## F-statistic: 119.4 on 1 and 82 DF,  p-value: < 2.2e-16

Para ajustarse R^2 , primero transformamos la variable predicha.

library(purrr)

y_predicted <- model_exponential$fitted.values %>% 
  map_dbl(~exp(.+0.09939^2/2))

y_predicted
##        1        2        3        4        5        6        7        8 
## 4.268591 4.289894 4.309227 4.330733 4.351647 4.373365 4.394485 4.416416 
##        9       10       11       12       13       14       15       16 
## 4.438457 4.459892 4.482149 4.503795 4.526272 4.548861 4.569361 4.592165 
##       17       18       19       20       21       22       23       24 
## 4.614342 4.637370 4.659766 4.683021 4.706392 4.729121 4.752722 4.775674 
##       25       26       27       28       29       30       31       32 
## 4.799508 4.823461 4.845198 4.869379 4.892894 4.917313 4.941060 4.965719 
##       33       34       35       36       37       38       39       40 
## 4.990502 5.014602 5.039628 5.063966 5.089238 5.114637 5.138512 5.164156 
##       41       42       43       44       45       46       47       48 
## 5.189096 5.214993 5.240177 5.266329 5.292612 5.318171 5.344712 5.370523 
##       49       50       51       52       53       54       55       56 
## 5.397326 5.424262 5.448707 5.475899 5.502344 5.529804 5.556509 5.584240 
##       57       58       59       60       61       62       63       64 
## 5.612109 5.639211 5.667355 5.694724 5.723144 5.751707 5.777627 5.806462 
##       65       66       67       68       69       70       71       72 
## 5.834503 5.863621 5.891938 5.921342 5.950894 5.979632 6.009474 6.038496 
##       73       74       75       76       77       78       79       80 
## 6.068632 6.098918 6.126404 6.156978 6.186712 6.217588 6.247614 6.278794 
##       81       82       83       84 
## 6.310129 6.340603 6.372247 6.403020

Y luego ajustamos R^2:

r_squared_exponential <- (cor(y_predicted,gasoline_df$gasoline))^2

adj_r_squared_exponential <- 1-(1-r_squared_exponential)*
  ((nrow(gasoline_df)-1)/(nrow(gasoline_df)-1-1))
adj_r_squared_exponential
## [1] 0.6476753

A veces, una serie temporal cambia la dirección, este tipo de variable objetivo se ajusta a los modelos de tendencia polinomiales . Si hay un cambio de dirección, usamos el modelo de tendencia cuadrática:

y_{t}={0}+{1}t+_{2}t^2+

Si hay dos cambios de dirección, usamos los modelos de tendencia cúbica:

y_{t}={0}+{1}t+{2}t^2+{3}t^3+

Estacionalidad

La estacionalidad representa las repeticiones en un período específico de tiempo. Las series temporales con observaciones semanales, mensuales o trimestrales tienden a mostrar variaciones estacionales que se repiten cada año.

Análisis de descomposición

son la tendencia, estacionalidad y los componentes aleatorios de la serie. Cuando la variación estacional aumenta a medida que aumenta la serie temporal, usaremos el modelo multiplicativo.

#we create a time series variable for the plot
gasoline_ts <- ts(gasoline_df$gasoline,start = c(2013,1),end = c(2019,12),frequency = 12)
#The plot have all the components of the series(T, S, I)
plot(gasoline_ts,lwd=2,ylab="Gasoline")

Como podemos ver en la gráfica anterior, parece que el modelo multiplicativo sería más adecuado para los datos. Cuando analizamos el gráfico, vemos algunas variaciones en la coyuntura causadas por la expansión o contracción de la economía. A diferencia de la estacionalidad, las fluctuaciones de la coyuntura pueden durar varios meses o años.

Extrayendo la estacionalidad

Los promedios móviles ( ma ) generalmente se usan para separar el efecto de una tendencia de la estacionalidad.Utilizamos el promedio móvil acumulativo ( CMA ), que es el promedio de dos promedios consecutivos, para mostrar el promedio móvil de orden par.

El valor predeterminado del argumento ‘centro’ de la función ma sigue siendo VERDADERO para obtener el valor CMA.

library(forecast)
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff
gasoline_trend <- forecast::ma(gasoline_ts,12)
#tambien se puede realizar de esta manera:
gasoline_detrend <- gasoline_ts/gasoline_trend

El promedio aritmético se usa para determinar el valor común para cada mes. Al hacer esto, eliminamos el componente aleatorio y restamos la estacionalidad de la variable de tendencia. Se llama índice estacional no ajustado:

unadjusted_seasonality <- sapply(1:12, function(x) mean(window(gasoline_detrend,c(2013,x),c(2019,12),deltat=1),na.rm = TRUE)) %>% round(4)

Los índices estacionales de datos mensuales deben completarse a 12, con un promedio de 1; Para hacerlo, cada índice estacional no ajustado se multiplica por 12 y se divide por la suma de 12 índices estacionales no ajustados:

adjusted_seasonality <- (unadjusted_seasonality*(12/sum(unadjusted_seasonality))) %>% 
                        round(4)

El promedio de todos los índices estacionales ajustados es 1; Si un índice es igual a 1, eso significa que no habría estacionalidad. Para trazar la estacionalidad realizamos lo siguiente:

#building a data frame to plot in ggplot
adjusted_seasonality_df <- data_frame(
  months=month.abb,
  index=adjusted_seasonality)
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
#converting char month names to factor sequentially
adjusted_seasonality_df$months <- factor(adjusted_seasonality_df$months,levels = month.abb)

ggplot(data = adjusted_seasonality_df,mapping = aes(x=months,y=index))+
 geom_point(aes(colour=months))+
 geom_hline(yintercept=1,linetype="dashed")+
 theme(legend.position = "none")+
 ggtitle("Seasonality of Unleaded Gasoline Prices for 95 Octane")+
 theme(plot.title = element_text(h=0.5))

Como podemos ver hay una disminución estacional del 3% en enero y diciembre, ademas hay un aumento estacional del 2% en mayo y junio. La estacionalidad no se ve en marzo, julio y agosto; porque sus valores de índice son aproximadamente iguales a 1.

Descomponiendo las series temporales gráficamente

Primero mostraremos la línea de tendencia en la serie temporal:

#Trend is shown by red line
plot(gasoline_ts,lwd=2,ylab="Gasoline")+
lines(gasoline_trend,col="red",lwd=3)

## integer(0)

Y ahora aislamos la línea de tendencia de la trama del grafico:

plot(gasoline_trend,lwd=3,col="red",ylab="Trend")

Observando la línea de estacionalidad:

plot(as.ts(rep(unadjusted_seasonality,12)),lwd=2,ylab="Seasonality",xlab="")

Por ultimo, mostraremos la línea de aleatoriedad:

randomness <- gasoline_ts/(gasoline_trend*unadjusted_seasonality)
plot(randomness,ylab="Randomness",lwd=2)

Agradecimientos

Autor: Selcuk Disci