El conjunto de datos que se va utilizar corresponde al precio de la gasolina en Turquía entre los años 2013 y 2020. El dataset se encuentra aquí.

library(readxl)

gasoline_df <- read_excel("gasoline_df.xlsx")

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

Primero queremos ver cómo fueron los precios en el periodo y analizar las líneas de tendencia ajustada. Hay que instalar y ejecutar la librería tidyverse. Como la variable de respuesta se mide logarítmicamente, tenemos que obtener el exponente de la misma y crear un dataframe diferente para la línea de tendencia exponencial.

library(tidyverse)

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

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 podemos observar en el gráfico los modelos cúbicos y exponenciales se superponen y parecen ser los que mejor se ajustan a los datos, pero vamos a analizar cada método en detalle.

Modelos predictivos de tendencia

Crearemos el modelo de tendencia lineal \(y_t=\beta_0 + \beta_1 t + \epsilon_t\). Para comparar los modelos, extraemos los coeficientes de determinación ajustados de cada modelo de tendencia.

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

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

adj_r_squared_linear
  [1] 0.6064

Ahora creamos el modelo exponencial \(ln(y_t)=\beta_0 + \beta_1 t + \epsilon_t\) (también ser podría haber utilizado exp.model anteriormente usado). Para saber cuál es error y el nivel de significación de los coeficientes realizamos un summary del modelo.

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

Como el p-valor de los coeficientes y del modelo son menos de 0.05, son significativos en un 5%. A continuación ajustaremos el \(R^2=1-(1-R^2)(\frac{n-1}{n-k-1})\) transformando la variable objetivo.

library(purrr)

y_predicted <- model_exponential$fitted.values %>% 
  map_dbl(~exp(.+0.09939^2/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 de dirección por muchas razones. El modelo que mejor ajusta esto son los modelos de tendencia polinómica. Si hay un cambio de dirección usamos el modelo de tendencia cuadrática \(y_t=\beta_0 + \beta_1 t + \beta_2 t^2 + \epsilon\). Si hay dos cambios de dirección usamos el modelo de tendencia cúbico \(y_t=\beta_0 + \beta_1 t + \beta_2 t^2 + \beta_3 t^3 + \epsilon\). Ahora, ejecutamos el mismo proceso que utilizamos con el modelo lineal.

model_quadratic <- lm(data = gasoline_df,gasoline_df$gasoline~poly(gasoline_df$date,degree=2))
model_cubic <- lm(data = gasoline_df,gasoline_df$gasoline~poly(gasoline_df$date,degree=3))

adj_r_squared_quadratic <- summary(model_quadratic) %>% 
  .$adj.r.squared

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

adj_r_squared_quadratic;adj_r_squared_cubic
  [1] 0.8659978
  [1] 0.8652594

Para comparar modelos, creamos una variable que compare los \(R^2\).

adj_r_squared_all <- c(linear=round(adj_r_squared_linear,4),
                   exponential=round(adj_r_squared_exponential,4),
                   quadratic=round(adj_r_squared_quadratic,4),
                   cubic=round(adj_r_squared_cubic,4))

adj_r_squared_all
       linear exponential   quadratic       cubic 
       0.6064      0.6477      0.8660      0.8653

Como podemos observar, los resultados justifican el gráfico que hemos realizado más arriba, los modelos polinómicos son mejores que los modelos lineales y los modelos exponenciales. A pesar de que los modelos polinómicos paracen iguales, el modelo cuadrático es un poco mejor que el modelo cúbico.

Estacionalidad

El componente de estacionalidad representa las repeticiones en un periodo específico de tiempo.

Análisis de descomposición: \(T\), \(S\), \(I\) son tendencia, estacionalidad y componentes aleatorios de la serie respectivamente. Cuando la variación estacional aumenta al mismo tiempo que la serie temporal aumenta, usamos el modelo multiplicativo \(y_t=T_t \times S_t \times I_t\). Si la variación parece constante, usaremos los modelos aditivos \(y_t=T_t + S_t + I_t\). Para saber qué modelo ajusta, tenemos que mirar en el gráfico.

gasoline_ts <- ts(gasoline_df$gasoline,start = c(2013,1),end = c(2019,12),frequency = 12)

plot(gasoline_ts,lwd=2,ylab="Gasoline")

Como podemos observar en el gráfico de arriba, parece que el modelo multiplicativo se ajusta mejor a los datos.

Extrayendo la estacionalidad

Las medias móviles se suelen usar para separar el efecto de la tendencia de la estacionalidad \(Media \hspace{2mm} movil=\frac{Media \hspace{2mm} de \hspace{2mm} las \hspace{2mm} ultimas \hspace{2mm} m \hspace{2mm} observaciones}{m}\)

library(forecast)

gasoline_trend <- ma(gasoline_ts, order=12)

gasoline_detrend <- gasoline_ts/gasoline_trend

gasoline_detrend
             Jan       Feb       Mar       Apr       May       Jun       Jul
  2013        NA        NA        NA        NA        NA        NA 1.0159944
  2014 1.0088095 1.0202510 1.0208072 1.0202957 1.0160643 1.0205456 1.0392090
  2015 0.8945041 0.9449850 0.9724972 0.9930729 1.0458051 1.0775966 1.0522201
  2016 0.9685894 0.9412665 0.9758898 0.9915278 1.0220963 1.0035468 0.9584137
  2017 1.0570470 1.0593789 1.0153846 1.0231804 0.9854050 0.9664046 0.9601515
  2018 0.9921295 0.9683846 0.9865800 0.9881865 1.0069930 1.0064516 1.0006682
  2019 0.9032014 0.9461176 1.0008177 1.0486684 1.0567367 1.0183727        NA
             Aug       Sep       Oct       Nov       Dec
  2013 1.0130612 1.0039594 0.9887179 0.9750636 0.9891277
  2014 1.0331057 1.0354424 1.0216833 0.9834836 0.9327708
  2015 1.0086665 1.0124404 0.9951220 0.9800527 0.9757943
  2016 0.9391023 0.9848217 1.0024314 0.9872311 1.0341844
  2017 0.9761307 0.9860451 0.9942910 1.0156797 0.9933388
  2018 1.0512344 1.0866794 1.0787217 1.0092450 0.9240756
  2019        NA        NA        NA        NA        NA

Cada mes tiene múltiples ratios, donde cada ratio coincide con un año diferente. La media aritmética es usada para determinar el valor común para cada mes. Haciendo esto, eliminamos el componente aleatorio y sustraemos la estacionalidad de la tendencia. Esto se le 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)

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("Seasonality of Unleaded Gasoline Prices for 95 Octane")+
 theme(plot.title = element_text(h=0.5))

Como podemos observar, hay aproximadamente un 3% de disminución estacional en enero y diciembre y un 2% de incremento estacional en Mayo y Junio. No se observa estacionalidad en los meses de marzo, julio y agosto porque el valor de sus índices es cercano a 1.

Descomponiendo las series temporales gráficamente

Primero observamos la tendencia en la serie temporal.

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

  integer(0)

Y ahora solo la tendencia.

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

Vamos a ver la estacionalidad.

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

Y también observamos la aletoriedad \(I_t=\frac{y_t}{S_t \times T_t}\).

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


Referencias

Sanjiv Jaggia, Alison Kelly (2013). Business Intelligence: Communicating with Numbers.

Selcuk Disci (2020). Trend Forecasting Models and Seasonality with Time Series. https://www.r-bloggers.com/trend-forecasting-models-and-seasonality-with-time-series/