En este documento utilizaremos los datos de los precios en Turquía, ya que allí los precios de la gasolina son un problema y siempre hay debate en ellos. Utilizaremos los datos de la empresa gasolina 95 sin plomo de la empresa Petrol Ofisi, los datos estarán organizados mensualmente desde 2013 a 2017 (analizaremos el transcurso de estos 7 años) con la moneda turca (Lira).

#setwd("C:/Users/alber/OneDrive/Escritorio/Universidad/Asignaturas Universidad/3º año/Segundo cuatrimestre/Análisis estadístico de series económicas/Tema 3")
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
# fecha gasolina 
# 1 2013-01-01 4.67 
# 2 2013-02-01 4.85 
# 3 2013-03-01 4.75 
# 4 2013-04-01 4.61 
# 5 2013-05-01 4.64 
# 6 2013 -06-01 4.72

En primer lugar, veremos como van los índices de los precios de la gasolina durante los periodos y analizaremos las líneas de tendencias ajustadas.

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()
library(ggplot2)

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

#Líneas de tendencia

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)

Primero tendremos que crear el modelo exponencial ya que si no al representar no nos irá, ya que no tendremos el modelo exponencial. Para hacer el modelo exponencial tendremos que hacerlo con el logarítmico, donde tendremos que obtener el exponente de la misma y crear un marco de datos diferente para la línea de tendencia exponencial.

Como vemos en el gráfico anterior, el modelo cúbico y exponencial se superponen a la línea y por tanto se ajustan mejor a los datos. De todas formas, analizaremos cada modelo por separado para verlo mejor y en detalle.

\(-Modelos\) \(de\) \(tendencia\) \(de\) \(pronóstico:\)

\(-Tendencia\) \(lineal:\) \(y_t\), el valor de la serie en un momento dado,que sería t,donde se describe como: \(y_t\)=\(\beta_0\)+\(\beta_1\)t+\(\epsilon\), donde \(\beta_1\) y \(\beta_2\) son los coeficientes. En este caso, en R sería:

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

Para comparar los diversos modelos tendremos que extraer los coeficientes de determinación ajustados.

\(-Modelo equilibrado:\)

El modelo equilibrado en este caso será: \(R ^ 2\) = 1- (1-\(R ^ 2\)) ({n-1}\(\frac\){nk-1}) donde n será el tamaño de la muestra, k el número de variables y \(R^2\) será el coeficiente de determinación que este es el cuadrado de correlación de los valores de de observación con los valores predichos: \(r^2_y\hat{y}\)

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

\(-Tendencia\) \(Exponencial:\)

A diferencia de la tendencia lineal, permite que la serie aumente a un ritmo creciente en cada período, por tanto, el modelo será:

ln(\(y_t\))=\(\beta_0\)+\(\beta_1\)t+\(\epsilon_t\).

Para hacer predicciones sobre este modelo ajustado, utilizaremos la función exponencial:

\(\hat{y}_t\) = exp(\(b_ {0}\) + \(b_ {1}\) t + \(s_ {e}\) \(^ {2/2}\))

Para encontrar \(s_{e}\), ejecutamos la función como se muestra a continuación. Además, podemos ver el nivel de significación de los coeficientes y el modelo. Como vemos a continuación, debido a que el valor p de los coeficientes y el modelo son menores a 0.05, son significativos al nivel de significancia del 5%.

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 ajustar \(R^2\), primero transformaremos la variable predicha:

library (purrr) 

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

Y luego ejecutaremos la fórmula anterior que hemos puesto anteriormente para ajustar \(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))

A veces una serie temporal cambia de 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}\) = \(\beta_ {0}\) +\(\ beta_ {1}\)t + \(\beta_ {2}\)\(t^ 2\) + \(\epsilon\)

Si hay 2 cambias de dirección, usamos el modelo de tendencia cúbica:

\(y_ {t}\) = \(\beta_ {0}\) +\(\ beta_ {1}\)t + \(\beta_ {2}\)\(t^ 2\) +\(\beta_ {2}\)\(t^3\)+ \(\epsilon\)

Por tanto, ejecutaremos el mismo proceso que en proceso lineal:

#Model variables
#Model variables
model_quadratic <- lm(data =gasoline_df,gasoline~poly(date,2))
model_cubic <- lm(data = gasoline_df,gasoline~poly(date,3))

#adjusted coefficients of determination
adj_r_squared_quadratic <- summary(model_quadratic) %>% 
  .$adj.r.squared

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

Para comparar todos los modelos, creamos una variable que reune todas las variables \(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))

Como podemos ver a continuación, los resultados justifican el gráfico que creamos antes; Los modelos polinómicos son mucho mejores que los modelos lineales y exponenciales. A pesar de que los modelos polinómicos se ven casi iguales, el modelo cuadrático es ligeramente mejor que el modelo cúbico.

adj_r_squared_all
##      linear exponential   quadratic       cubic 
##      0.6064      0.6477      0.8660      0.8653

\(-Estacionalidad:\)

La estacionalidad representa las repeticiones en un periodo específico de tiempo. Las series de tiempo con observaciones semanales mensuales o trimestrales tienden a mostrar variaciones estacionales que se repiten cada año. Por ejemplo, la venta de artículos al por menor aumenta cada año en el período navideño o los recorridos festivos aumentan en el verano. Para detectar la estacionalidad, utilizaremos el análisis de descomposición.

\(-Análisis\) \(de\) \(descomposición:\)

T (tendencia), S (estacionalidad) y I (componentes aleatorias de la serie). Cuando la variación estacional aumenta a medida que aumenta la serie de tiempo, usaríamos el modelo multiplicativo.

\(y_t\)=\(T_t\)x\(S_t\)x\(I_t\)

Si la variación parece constante, deberíamos usar el modelo aditivo.

\(y_t\)=\(T_t\)+\(S_t\)+\(I_t\)

Para ver que el modelo se ajusta, tendremos que mirarlo en el gráfico:

#Crearemos una variable de serie temporal para el gráfico
gasoline_ts <- ts(gasoline_df$gasoline,start = c(2013,1),end = c(2019,12),frequency = 12)
#El gráfico tiene todos los componentes del diagrama (T, S, I)
plot(gasoline_ts,lwd=2,ylab="Gasoline")

Como podemos ver en el gráfico, el modelo multiplicativo es el que mejor se ajusta a los datos, especialmente si miramos del año 2016 al 2019. Analizando el gráfico, vemos alguna variación en la conyuntura causada por la expansión o contración de la economía. A diferencia de la estacionalidad, las fluctuaciones de la conyuntura puede durar varios meses o años.

Además. la magnitud de las fluctuaciones de la conyuntura pueden variar con el tiempo porque las fluctuaciones difieren en magnitud y longitud, por lo que son difíciles de reflejar con datos históricos.

\(-Extrayendo\) \(la\) \(estacionalidad:\)

Los promedios móviles (ma) generalmente se usan para separar el efecto de una tendencia de la estacionalidad.

m term moving averages=average of the last m observations\(/\)m.

Utilizaremos el promedio móvil acumulativo (CMA), que es el promedio de dos promedios consecutivos, que sirve 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)

\(\bar y_t\) elimina las variaciones estacionales y aleatorias, por lo tanto: \(\bar y_t\)=\(T_t\).

Ya que \(y_t\)=\(T_t\)x\(S_t\)x\(I_t\) e \(\bar y_t\)=\(T_t\) donde la variable tendencia se encuentra por \(\bar y_t\)\(/\)\(y_t\)=\(S_t\)x\(I_t\); que también llamaremos método de relación de promedios móviles.

gasoline_detrend <- gasoline_ts/gasoline_trend

Cada mes tiene múltiples razones, donde cada razón coincide con un año diferente. En esta muestra, cada mes tiene siete proporciones diferentes. 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, esta 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:

#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)
library(ggplot2)

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 en el gráfico anterior, hay aproximadamente una disminución estacional del 3% en enero y diciembre; por otro lado, 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\) \(de\) \(tiempo\) \(gráficamente\)

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

#La tendencia se muestra en la línea roja
plot(gasoline_ts,lwd=2,ylab="Gasoline")+
lines(gasoline_trend,col="red",lwd=3)

## integer(0)

Y separaremos la línea de tendencia de la trama:

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

Ahora veamos la línea de estacionalidad:

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

Y finalmente veremos la línea de aleatoriedad, que será la función \(I_t\)=\(y_t\)\(/\)(\(S_t\)x\(T_t\))

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

\(-Agradecimientos:\)

En este apartado queremos dar las gracias a Selcuk Disci por la publicación del artículo Trend Forecasting Models and Seasonality with Time Series en RPubs.