Este Trabajo es una adaptación de un artículo publicado por Selcuk Disci publicado en R-bloggers cuyo enlace encontraréis al final del artículo.

A los habitantes de Turquía les gusta mucho conducir, siempre cogen el coche 🚗 para ir de un lado al otro y por eso los precios de la gasolina ⛽ allí siempre es un gran dilema y los turcos siempre están quejandose de sus precios. En el estudio posterior se han utilizado datos de una compañía distribuidora de gasolina llamada Petrol Ofisi ,concretamente se han recogido los precios de gasolina 95 sin plomo entre los años 2013 y 2020 ,los datos fueron recogidos mensualmente y los precios son en Lira, moneda que se utiliza en Turquía. Puedes descargar los datos aqui

library(kableExtra)
kable(head(gasoline_df)) %>%
  kable_styling(bootstrap_options = c("condensed","bordered"),full_width = F, font_size = 14,position = "float_right")%>%
                 column_spec(column=1, bold = TRUE,  color = "black", background = "lightgrey")%>%
                   column_spec(column=2, bold = TRUE,  color = "black", background = "lightgrey")%>%
  row_spec(1:6,hline_after =T)
date gasoline
2013-01-01 4.67
2013-02-01 4.85
2013-03-01 4.75
2013-04-01 4.61
2013-05-01 4.64
2013-06-01 4.72



En la tabla se observan como han sido recogidos los datos. En la columna de la izquierda se representa el tiempo y en la de la derecha los precios de gasolina en la moneda turca.Antes de todo vamos a visualizar los datos y diferentes líneas de tendencias representadas sobre dichos datos. En el eje X representaremos las fechas y en eje Y los precios de la gasolina. Se intentará estudiar si hay tendemcia ,estacionalidad,ciclo en los precios de gasolina en Turquía.



Visualizar los datos y diferentes líneas de tendencias


Vamos a representar los datos y ver si hay alguna línea (tendencia) que se ajuste a estos datos y para ello representaremos múltiples tipos de modelos como pueden ser:lineales,exponenciales…

#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)))
#Trend lines 
library(tidyverse)
library(emoGG)
grafico=ggplot(gasoline_df, aes(x = date, y = gasoline))+add_emoji(emoji="26fd") + geom_point(color = "black") +
  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) 

grafico + theme(
  plot.background = element_rect(fill = "#F3F6FA"),
  panel.background = element_rect(fill = '#E8E8E8',
                                colour = '#C0C0C0'
                                ),
  legend.background = element_rect(fill = '#F3F6FA'),
  legend.key = element_rect(fill="#F3F6FA", color="#F3F6FA"),
                                  
  panel.grid.major = element_line(size = 0.5, linetype = 'solid',
                                colour = "#CECACA"), 
  panel.grid.minor = element_line(size = 0.25, linetype = 'solid',
                                colour = "#CECACA")
)

Se puede observar que el modelo cúbico y exponencial se solapan. Si observamos el gráfico de arriba,son los que mejor se están ajustando. Pero no nos vamos a limitar a ese gráfico sino que haremos un análisis a continuación de los distintos modelos.


Modelos de Tendencia

Tendencia Lineal : \(\color{#69FFA4}{y_{t}=\beta_{0}+\beta_{1}t+\epsilon_{t}}\)   , \(\color{#91FFEB}{\beta_{0}}\) y \(\color{#91FFEB}{\beta_{1}}\) son coeficientes   y   \(\color{#FFFF6E}\epsilon\)   el error.

Vamos a crear una variable donde vamos a guardar nuestro modelo lineal. A posteriori sacaremos el coeficiente de determinación ajustado, que nos ayudará a comparar diferentes modelos y nos guiará hacía el mejor modelo.

modelo_linear <- lm(data = gasoline_df,gasoline~date)

A continuación se recuerda la fórmula que se emplea para sacar coeficiente de determinación ajustado:

Coeficiente Ajustado: \(\color{#FFA0B9}{R^2=1-(1-R^2)(\frac{n-1}{n-k-1})}\)
\(\color{#FFEFA0}n:\)Tamaño Muestral
\(\color{#FF91C8}k:\) Número de Variables
\(\color{#69FFA4}{R^2}:\) coeficiente de determinación que es corelación al cuadrado de los valores observados con los valores predichos. \(\color{#69FFA4}{r^2_{y\hat{y}}}\)

adj_r_cuadrado_linear <- summary(modelo_linear) %>% 
  .$adj.r.squared %>% 
  round(4)

%>% es un operador de r que nos permite concatenar varias operaciones, para más información os dejo el Enlace
Tendencia Exponencial: a diferencia de la tendencia lineal, el método exponencial permite que la serie aumente a un ritmo creciente en cada período,y se describe como: \(\color{#FAC5DE}{\ln(y_{t})=\beta_{0}+\beta_{1}t+\epsilon_{t}}\).

Para realizar predicciones utilizando un modelo ajustado primero debemos transformar el modelo debido a que los datos originales fueron transformados mediante logaritmo natural \(\color{#2BFFFC}{\ln(y_{t})}\) y para invertir \(\color{#A4DDF5}\ln\) utilizamos \(\color{#A4DDF5}{exp}\).  Entonces, el modelo quedaría de la siguiente forma: \(\color{#FCD6FF}{\hat{y_{t}}=exp(b_{0}+b_{1}t+s_{e}^2/2)}\)  sumamos la mitad de desviación estándar al cuadrado para que \(\color{#EEF5A4}{\hat{y_{t}}}\) no sea subestimada por \(\color{#EEF5A4}{y_{t}}\). para extraer \(\color{#EEF5A4}{s_{e}}\) utilizamos el commando summary(modelo exponencial). Para ello primero vamos a definir modelo exponencial.

modelo_exponential <- lm(data = gasoline_df,log(gasoline)~date)

summary(modelo_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.0758581  0.2495848  -4.311  4.5e-05 ***
## date         0.0001606  0.0000147  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

Podemos observar en el summary que todos los coeficientes son significativos a un nivel de 5% porque el pvalor es menor a <0.05 de todos los coeficientes.

Para Obtener \(\color{#FFC836}{R^2}\) ajustado primero debemos transformar nuestros valores predichos por el modelo exponencial y después con la formula presentada arriba se calcula \(\color{#FFC836}{R^2}\) ajustado.

library(purrr)

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

Y ahora con la formula se calcula \(\color{#FFC836}{R^2}\) ajustado.

r_cuadrado_exponential <- (cor(y_predichos,gasoline_df$gasoline))^2

adj_r_cuadrado_exponential <- 1-(1-r_cuadrado_exponential)*
  ((nrow(gasoline_df)-1)/(nrow(gasoline_df)-1-1))

A veces una serie temporal cambia de dirección, es decir, sube y luego baja de manera alternativa. Dicho cambio puede ser por varias razones y es por eso que este tipo de serie la trataremos por los modelos polinómicos.

Si solo hay un cambio en la dirección se utilizará el modelo de Tendencia cuadrático : \(\color{#2BFFFC}{y_{t}=\beta_{0}+\beta_{1}t+\beta_2{t}^2+\epsilon_{t}}\)

y si hay 2 cambios en la dirección utilizaremos modelos de Tendencia Cúbicos : \(\color{#37FF7D}{y_{t}=\beta_{0}+\beta_{1}t+\beta_2{t}^2+\beta_3{t}^3+\epsilon_{t}}\)

Y a continuación se hará el análisis de igual manera que en los otros modelos, es decir, los modelos serán guardados en variables y luego se sacará el coeficiente de determinación \(\color{#FFC836}{R^2}\) ajustado para comparar los modelos.

#variables donde se guardan los modelos
model_cuadratic <- lm(data = gasoline_df,gasoline~poly(date,2))
model_cubic <- lm(data = gasoline_df,gasoline~poly(date,3))
#coeficientes de determinación ajustados 
adj_r_cuadrado_cuadratic <- summary(model_cuadratic) %>% 
  .$adj.r.squared

adj_r_cuadrado_cubic <- summary(model_cubic) %>% 
  .$adj.r.squared

Ahora, para comprobar los coeficientes de determinación ajustados de todos los modelos ajustados, vamos a crear una tabla con \(\color{#FFC836}{R^2}\) el valor ajustado de cada modelo.

adj_r_cuadrado_todo <- data.frame(y=c('Lineal','exponencial'
                             ,'cuadrático','cúbico'),
                             x=c(round(adj_r_cuadrado_linear,4),
                             round(adj_r_cuadrado_exponential,4),
                             round(adj_r_cuadrado_cuadratic,4),
                             round(adj_r_cuadrado_cubic,4)))
kable(adj_r_cuadrado_todo,col.names = c('Modelo','$R^2$')) %>%
  kable_styling(bootstrap_options = c("condensed","bordered"),full_width = F, font_size = 14,position = "float_left")%>%
                 column_spec(column=1:2, bold = TRUE,  color = "black", background = "lightgrey")
Modelo \(R^2\)
Lineal 0.6064
exponencial 0.6477
cuadrático 0.8660
cúbico 0.8653


Como se puede observar, los modelos Polinómicos son los que mejor se ajustan a estos datos y los modelos lineales y exponencial no se ajustan tan bien. Esta conclusión también se veía reflejada en el gráfico y par dar más validez a los argumentos,vemos que los coeficientes de determinación ajustados apuntan hacía la misma dirección. Entre modelo cuadrático y el cúbico no hay mucha diferencia, aunque el cuadrático es un poco mejor.


Estacionalidad

Se llama estacionalidad a los patrones de comportamiento que regularmente exhibe una variable en momentos específicos del año. Por ende, la conducta de estas variables se convierte en predecible y al conocer esta característica, su estudio se facilita. Un ejemplo de este concepto sería la evolución del consumo de abrigos en un año. Claramente, no podemos comparar el uso de estos en diciembre y en junio.

Descomposición de la serie : A una Serie temporal se le pueden extraer varias componentes entre las cuáles están:

\(T :\) Tendencia,  \(S :\) Estacionalidad,  \(I :\) Componente aleatoria(Error)

Si en una serie temporal la variación en estacionalidad va aumentando conforme se incrementa la serie, es decir, la estacionalidad se multiplica o es mayor en comparación con el período anterior, entonces en estos casos es preferible utilizar un modelo multiplicativo : \(\color{#26FF88} {y_t=T_t\times S_t\times I_t}\)

En cambio, si la variación en la estacionalidad se mantiene constante a lo largo de la serie, es preferible utilizar el modelo aditivo: \(\color{#26E1FF}{y_t=T_t+ S_t+ I_t}\)

Para saber que modelo se ajusta mejor debemos visualizar la serie y ver si hay indicios de estacionalidad y si esta es multiplicativa,incrementa a lo largo de la serie, o se mantiene constante.

#Vamos a crear una serie temporal para posterior representarla
gasoline_ts <- ts(gasoline_df$gasoline,start = c(2013,1),end =
                    c(2019,12),frequency = 12)
#Gráfico muestra todas las componentes(T, S, I)
plot(gasoline_ts,lwd=2,ylab="Gasoline")

#con ggplot
p <- ggplot(data = gasoline_df, aes(x = date, y = gasoline)) + 
     geom_line(color = "#00AFBB", size = 1.5)
p

Como podemos observar en los gráficos, el mejor modelo para analizar sería el modelo multiplicativo, debido a que después del año 2016 podemos observar que la serie incrementa, pero además hay que añadir que a esta variación se le añade que la variación anual también aumenta, no se mantiene constante.

También en el gráfico,podemos observar que la serie tiene una serie de subidas y bajadas,y algunas de esas bajadas (subidas) duran varios meses, llegando incluso a un año. Esto es debido a las variaciones en la economía que también se denominan ciclos.

Desestacionalizar

Medias móviles son las que se utilizan habitualmente para separar la tendencia de la componente estacional. Así las medias móviles son una lista de números en la cuál cada uno es el promedio de un subconjunto de datos de los originales.El método de medias móviles suaviza la serie de datos original, considerándose esta nueva serie suavizada la tendencia.

Por ejemplo, si se tiene un conjunto de 100 datos, el primer valor de la serie de medias móviles podría ser el promedio de los primeros 25 términos, luego el promedio de los términos 2 al 26, el tercer elemento de los términos 3 al 27 y así, hasta por último el promedio de los últimos 25 números del 76 al 100.

\(n\,termino\,media\,movil=\frac{media\,de \,los\,ultimos\,terminos}{n}\)

Utilizaremos las Medias Móviles acumuladas(CMA) que es la media de medias móviles.Es decir primero cálculamos medias móviles y luego haremos su media de dos a dos.

Utilizaremos la librería forecast para extraer las componentes de una serie temporal.

library(forecast)

gasoline_trend <- forecast::ma(gasoline_ts,12)

1.primero quitaremos la tendencia.

gasoline_sintendencia <- gasoline_ts/gasoline_trend

2.Segundo calculamos los índices de estacionalidad. Como la media móvil se calculan como media de n términos consecutivos ,entonces si en una serie temporal hay varios años, entonces a cada mes le corresponderán varias medias. Pero cada uno corresponderá a un año diferente. Por ejemplo, en nuestro caso cada mes tiene 7 medias.
Haremos una media aritmética para que cada mes tenga un único valor común.De esta manera se podrá representar y comparar con valores de otros meses y así estudiar la estacionalidad de cada mes y con eso también eliminamos la componenete aleatoria y quitamos estacionalidad a la variable, ya sin tendencia creada anteriormente. Se le conoce como índice de estacionalidad no ajustado(unadjusted seasonal index) o también se le conoce como índices de variación estacional.

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

Los índices de estacionalidad básicamente nos indican ese mes o cuatrimestre como ha actuado, comparando con la estacionalidad media del año entero. Como la estacionalidad media del año es 1, por eso lo comparamos los índices de cada mes. Si son distintos a uno, significa ha habido una estacionalidad. Por ejemplo, si para un mes el indice de estacionalidad es de 1.3, eso quiere decir que ese mes ha tenido un 30% más de valor que el respectivo valor medio. Puediendo ser ese valor la renta, beneficios,etc.

Como vamos a sacar índices de variacion estacional por meses, en total los índices deben sumar 12 porque hay 12 meses. Para ajustar nuestros índices, de manera que den 12, cada índice estacional no ajustado se multiplica por 12 y se divide por la suma de los 12 índices estacionales no ajustados.

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

Ahora que tenemos los índices, vamos a realizar un plot para observar en qué meses la estacionalidad aumenta o decrece.

#building a data frame to plot in ggplot
adjusted_seasonality_df <- data_frame(
  months=month.abb,
  index=adjusted_seasonality)

#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,size=1.0))+
 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))


En el gráfico se observa que la línea horizontal representa al número 1, que nos indica que los puntos (meses) que esten sobre esa línea no presentarán estacionalidad, que son los meses de julio,agosto y marzo. Mientras que en los meses de enero y deciembre hay una disminución de estacionalidad en un 3% y los meses de mayo, junio y septiembre hay un incremento de 2% en la estacionalidad. En palabras simples, cuando hay incremento de los precios de gasolina en turquía, suben en comparación con los precios medios de todo el año ,mientras que cuando hay disminución estos tienden a bajar.

Extrayendo Componentes de la Serie Gráficamente


Primero mostraremos la tendencia calculada anteriormente sobre el gráfico de los datos y luego mostraremos un gráfico sólo con la tendencia.

gasoline_dummy=data.frame(gasoline_trend)
gasoline_trend1=gasoline_df
gasoline_trend1[,2]=gasoline_dummy[,1]
ggplot() + 
  geom_line(data = gasoline_df, aes(x = date, y = gasoline), color = "#4A50FE",size=1.5) +
  geom_line(data = gasoline_trend1,aes(x= date,y= gasoline) ,color = "red",size=1.5) +
  xlab('fecha') +
  ylab('precios gasolina')

#Tendencia
ggplot() + 
  geom_line(data = gasoline_trend1,aes(x= date,y= gasoline) ,color = "red",size=1.5) +
  xlab('fecha') +
  ylab('precios gasolina')

#Estacionalidad
Estacionalidad=data.frame(as.ts(rep(unadjusted_seasonality,12)))
colnames(Estacionalidad)<-'estacionalidad'
ggplot() + 
  geom_line(data = Estacionalidad, aes(x=1:144,y = estacionalidad), 
            color = "#555555",size=1.5,)+
  xlab(' ')


Y para terminar vamos a representar la componente aleatoria: \(\color{#7CFF90}{I_{t}=\frac {y_{t}}{S_{t}\times T_{t}}}\)

aleatorio <- gasoline_ts/(gasoline_trend*unadjusted_seasonality)
gasoline_dummy=data.frame(aleatorio)
aleatorio=gasoline_df
aleatorio[,2]=gasoline_dummy[,1]
ggplot() + 
  geom_line(data = aleatorio,aes(x= date,y= gasoline) ,color = "#2E1C7E",size=1.5) +
  xlab('fecha') +
  ylab('precios gasolina')