El siguiente informe corresponde a los precios de gasolina de 95 octanos sin plomo (medido en liras) en el estado de Turquía. Los datos han sido recogidos a lo largo de siete años, desde el año 2013 hasta 2020.

Cargamos el dataframe

#install.packages("readxl")
library(readxl)
gasolina <- read_excel("gasoline_df.xlsx", col_types = c("text", "date", "numeric"))

Renombramos los títulos de las columnas, convertimos los datos de la columna fecha en formato “fecha” y visualizamos las primeras observaciones de la tabla.

names(gasolina)[names(gasolina) == "...1"] <- "Observaciones"
gasolina$date=as.Date(gasolina$date, format = "%d/ %m/ %Y")
names(gasolina)[names(gasolina) == "date"] <- "Fecha"
names(gasolina)[names(gasolina) == "gasoline"] <- "Precio"
head(gasolina)
## # A tibble: 6 x 3
##   Observaciones Fecha      Precio
##   <chr>         <date>      <dbl>
## 1 1             2013-01-01   4.67
## 2 2             2013-02-01   4.85
## 3 3             2013-03-01   4.75
## 4 4             2013-04-01   4.61
## 5 5             2013-05-01   4.64
## 6 6             2013-06-01   4.72

Visualización de líneas de tendencia ajustadas

Antes de comenzar el análsis en sí, extraemos la parte exponencial de los datos en la columna “Precio”. Por lo tanto, creamos nuevas líneas exponenciales de tendencia.

exp.model <- lm(log(Precio)~Fecha,data = gasolina) 
exp.model.df <- data.frame(x=gasolina$Fecha,
                           y=exp(fitted(exp.model)))

Visualizamos gráficamente qué modelo se ajusta mejor a los datos. En cualquier caso, se analizará caso a caso y determinará el resultado analíticamente.

#install.packages("tidyverse")
library(tidyverse)
ggplot(gasolina, aes(x = Fecha, y = Precio)) + geom_point() +
  stat_smooth(method = 'lm', aes(colour = 'Lineal'), se = FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,2), aes(colour = 'Cuadrática'), se= FALSE) +
  stat_smooth(method = 'lm', formula = y ~ poly(x,3), aes(colour = 'Cúbica'), se = FALSE)+
  stat_smooth(data=exp.model.df, method = 'loess',aes(x,y,colour = 'Exponencial'), se = FALSE)

Observamos gráficamente que los modelos cuadrático y cúbico son extremadamente similares y poseen mayor significación respecto al resto.

A continuación, para escoger el mejor modelo, compararemos sus coeficientes de determinación ajustados, indicador del porcentaje total de comportamiento de las variables explicado por el modelo.

Modelo lineal

Creamos el propio modelo lineal y extraemos el comparador que nos interesa, el R^2 ajustado.

model_linear <- lm(data = gasolina,Precio~Fecha)
summary(model_linear)
## 
## Call:
## lm(formula = Precio ~ Fecha, data = gasolina)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04367 -0.43706 -0.01212  0.47015  0.90432 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -9.941e+00  1.341e+00  -7.414    1e-10 ***
## Fecha        8.962e-04  7.895e-05  11.352   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.534 on 82 degrees of freedom
## Multiple R-squared:  0.6111, Adjusted R-squared:  0.6064 
## F-statistic: 128.9 on 1 and 82 DF,  p-value: < 2.2e-16
adj_r_squared_linear <- summary(model_linear) %>% 
  .$adj.r.squared %>% 
  round(4)
adj_r_squared_linear
## [1] 0.6064

Modelo exponencial

Ejecutamos los comandos respectivos al modelo exponencial y transformamos la ecuación para predecir mejor con el modelo ajustado y hallamos R^2 ajustado.

model_exponential <- lm(data = gasolina,log(Precio)~Fecha)
summary(model_exponential)
## 
## Call:
## lm(formula = log(Precio) ~ Fecha, data = gasolina)
## 
## 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 ***
## Fecha        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
#install.packages("purrr")
library(purrr)
y_predicted <- model_exponential$fitted.values %>% 
  map_dbl(~exp(.+0.09939^2/2))
r_squared_exponential <- (cor(y_predicted,gasolina$Precio))^2
adj_r_squared_exponential <- 1-(1-r_squared_exponential)*
  ((nrow(gasolina)-1)/(nrow(gasolina)-1-1))
adj_r_squared_exponential
## [1] 0.6476753

Modelos polinómicos

Observando los cambios de dirección en la serie temportal en el apartado ‘’Visualización de líneas de tendencia ajustadas’’, es conveniente crear modelos polinómico. El primero de ellos será cuadrático, de grado 2.

model_quadratic <- lm(data =gasolina,Precio~poly(Fecha,2))
summary(model_quadratic)
## 
## Call:
## lm(formula = Precio ~ poly(Fecha, 2), data = gasolina)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5792 -0.2311 -0.0135  0.2427  0.8070 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      5.26560    0.03399  154.90   <2e-16 ***
## poly(Fecha, 2)1  6.06138    0.31155   19.46   <2e-16 ***
## poly(Fecha, 2)2  3.93929    0.31155   12.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3116 on 81 degrees of freedom
## Multiple R-squared:  0.8692, Adjusted R-squared:  0.866 
## F-statistic: 269.2 on 2 and 81 DF,  p-value: < 2.2e-16
adj_r_squared_quadratic <- summary(model_quadratic) %>% 
  .$adj.r.squared
adj_r_squared_quadratic
## [1] 0.8659978

El siguiente modelo polinómico a crear será el cúbico, de grado 3.

model_cubic <- lm(data = gasolina,Precio~poly(Fecha,3))
summary(model_cubic)
## 
## Call:
## lm(formula = Precio ~ poly(Fecha, 3), data = gasolina)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56860 -0.22494 -0.03312  0.24529  0.78562 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      5.26560    0.03409 154.478   <2e-16 ***
## poly(Fecha, 3)1  6.06138    0.31241  19.402   <2e-16 ***
## poly(Fecha, 3)2  3.93929    0.31241  12.609   <2e-16 ***
## poly(Fecha, 3)3 -0.23297    0.31241  -0.746    0.458    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3124 on 80 degrees of freedom
## Multiple R-squared:  0.8701, Adjusted R-squared:  0.8653 
## F-statistic: 178.7 on 3 and 80 DF,  p-value: < 2.2e-16
adj_r_squared_cubic <- summary(model_cubic) %>% 
  .$adj.r.squared
adj_r_squared_cubic
## [1] 0.8652594

Comparación de coeficientes de determinación ajustados

Antes de comparar, observamos que los p-valores de los cuatro modelos son inferiores a 0.05; por lo tanto, todos los modelos son significativos a un nivel de significación del 5% y aptos para la comparación intermodelos.

adj_r_squared_all <- c(lineal=round(adj_r_squared_linear,4),
                       exponencial=round(adj_r_squared_exponential,4),
                       cuadratico=round(adj_r_squared_quadratic,4),
                       cubico=round(adj_r_squared_cubic,4))
adj_r_squared_all
##      lineal exponencial  cuadratico      cubico 
##      0.6064      0.6477      0.8660      0.8653

Concluimos las mismas hipótesis que observamos en el análisis gráfico: los modelos polinómicos están por encima del resto y la diferencia que presentan entre ellos es mínima. Aún así, el modelo cuadrático predice un poco mejor que el cúbico.

Estacionalidad

La componente estacional de una variable es aquella que se repite cada x tiempo. Para trabajarla, utilizaremos análisis de descomposición del modelo. Si la estacionalidad aumenta según lo hace la serie, nos encontramos frente a un modelo multiplicativo. Por otra parte, si es constante, será un modelo aditivo.

Para averiguarlo, realizaremos un gráfico del precio de la gasolina en función del tiempo y lo comprobaremos.

gasoline_ts <- ts(gasolina$Precio,start = c(2013,1),end = c(2019,12),frequency = 12)
plot(gasoline_ts,lwd=2,ylab="Gasolina")

Como observamos, la estacionalidad no es la misma en todos los años; existe, pero en distintas magnitudes, por lo que utilizaremos un método multiplicativo.

Para extraer la estacionalidad, utilizaremos el método de medias móviles, que extrae el efecto de la tendencia del componente estacional.

#install.packages("forecast")
library(forecast)

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

Una vez controlada la variable tendencia, nos encontramos con los datos dependientes de estacionalidad y aleatoriedad. Recordemos la función del método multiplicativo:

.

Entonces, para extraer el componente aleatorio haremos el índice de estacionalidad sin ajustar; calcularemos la media artimética de los siete (número de años en la tabla) ratios de cada mes, para, de esta forma, eliminar la aleatoriedad y quedarnos con lo que deseamos, la estacionalidad.

Debido al planteamiento teórico de los índices, éstos se multiplicarán por 12 y dividirán por la suma de los doce índices estacionales no ajustados del año para que la media total sea 1.

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

adjusted_seasonality_df <- data_frame(
  months=month.abb,
  index=adjusted_seasonality)

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("Estacionalidad de los precios de la Gasolina 95 octanos sin plomo")+
  theme(plot.title = element_text(h=0.5))

Observamos cómo la estacionalidad es inferior a 1, hasta un 3%, en los meses invernales. En otras palabras, el precio se ve reducido los meses de invierno en el hemisferio norte. Por su parte, aumenta en mayo, junio y septiembre, los meses que rodean julio y agosto, conocidos por su gran capacidad turística en las tierras turcas. Otro dato relevante es la falta de estacionalidad en los propios meses mencionados anteriormente.

Finalizando, mostraremos la línea de tendencia junto a la serie para posteriormente enseñarla sola.

plot(gasoline_ts,lwd=2,ylab="Gasolina")+
  lines(gasoline_trend,col="red",lwd=3)

## integer(0)
plot(gasoline_trend,lwd=3,col="red",ylab="Tendencia")

Por último, los dos últimos gráficos representan la línea de estacionalidad y aleatoriedad, respectivamente. El primer gráfico es una serie del observado en los pasos previos.

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

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

Agradecimientos

Este trabajo ha podido ser realizado gracias al material compartido por el profesor Xavier Barber de la Universidad Miguel Hernández de Elche (España) y el informe previo realizado por Selcuk Disci, publicado en R-Bloggers el 2 de marzo de 2020. Puedes consultar su trabajo aquí.