Hands-on 2 ST_Regalos

Llamamos la siguientes librerias

library(ggplot2)
library(ggthemes)
library(dplyr)
library(lubridate)
library(forecast)
library(tidyverse)
library(janitor)

Datos del numero de veces que se busca la palabra regalos por semana en google

Componentes de la serie de tiempo

Cargamos datos

datos_regalos <- read.csv("data/datos_regalos.csv", skip = 1) %>% 
  janitor::clean_names()

La libreria janitor limpia los nombres Queremos transformar los datos a formato de serie de tiempo Usamos la funcion ts(), el argumento frequency en este caso semanal = 52 En el argumento start ponemos cuando obtubimos los datos 09/10/2016, semana=40

regalos.ts = ts(datos_regalos$regalos_mexico, start=c(2016,40), frequency = 52)
plot(regalos.ts,xlab="Tiempo",ylab="Regalos",main="N de veces que se busca Regalo")

Podemos observar

Varianza no constante Media no constante Estacionalidad,los picos ocurren en el mismo mes, el ciclo parece repetirse cada año.Ademas, la amplitud de las fluctuaciones aumenta en determinados meses Conlcuimos que la serie no es estacionaria.

Descomposicion de una serie en: Serie observada = Tendencia + Efecto estacional + Error.

regalos.ts.desc = decompose(regalos.ts)

Graficamos

plot(regalos.ts.desc, xlab='regalos')

Analicemos si el error que queda en la descomposicion tiene distribucion normal

Ho: tienen distribucion Normal (p value alfa=.05)

Ha: no tienen distribucion Normal (p value < alfa=.05)

qqline(y=as.integer(regalos.ts.desc$random))
shapiro.test(x=as.integer(regalos.ts.desc$random))

como el p value es muy pequeño concluimos que no tienen distribucion normal N(o,1)

ACF, PACF

Como es la autocorrelacion

Grafiquemos el correlograma de ACF, PACF

Si el coeficiente de autocorrelacion sobrepasa el limite entonces es significante

ACF

acf(regalos.ts)

Observamos 11 observaciones significantes

PACF

pacf(regalos.ts)

Observamos 4 observaciones significantes

Transformaciones

Aplicamos las transformaciones para eliminar la tendencia

Hacer la media contante y la varianza constante

Esto sera dificil debido a la estacionalidad,unos picos seran mes grandes

Diferencias de Box-Jenkins

Una forma sencilla de eliminar una tendencia aproximadamente lineal es diferenciar la serie, es decir, considerar la serie de diferencias entre una bservacion y la anterior en lugar de la serie original. Diferenciacion de series de tiempo lag con el que series de tiempo

diff.datos_regalos=diff(datos_regalos$regalos_mexico, lag=1)
plot(diff(datos_regalos$regalos_mexico, lag=1), type="l")

Desplazamos la serie

lag(datos_regalos$regalos_mexico, 1)

  [1]  NA  29  26  22  28  32  30  34  47  60  70  87  59  39  31  33  43
 [18]  53  72  75  22  23  21  20  22  25  26  22  23  33  44  51  29  25
 [35]  27  29  44  33  27  31  32  32  27  29  23  26  22  24  25  20  21
 [52]  28  29  30  32  25  27  30  37  40  46  64  73  87  57  39  35  36
 [69]  43  46  59  88  25  26  21  26  22  20  28  26  29  32  42  57  33
 [86]  25  32  28  43  26  25  28  28  26  27  26  24  24  23  28  24  24
[103]  24  28  32  30  24  27  27  31  32  43  51  75  87 100  71  42  36
[120]  36  41  42  53  99  24  22  21  21  20  24  29  24  25  30  37  60
[137]  31  23  24  31  45  33  27  30  25  20  24  25  21  25  22  22  23
[154]  22  21  24  25  26  27  25  23  31  35  37  50  59  74  87  76  39
[171]  35  33  36  42  49  97  26  20  21  21  15  13  15  17  18  18  32
[188]  57  34  24  22  28  33  51  31  24  25  28  24  25  25  19  24  22
[205]  23  20  22  19  22  23  25  24  24  24  31  38  41  51  60  83  99
[222]  36  40  34  32  36  43  62  35  20  20  18  21  19  20  19  23  30
[239]  35  50  37  21  21  23  27  46  27  20  25  25  22  20  20  19  22
[256]  22  20  20  17  18

Observamos que hemos eliminados la tencencia, haciendo media=0, constante datos_t= datos transformados

datos_t=data.frame(datos_regalos$semana[-1],diff.datos_regalos)

Cambiamos encabezados

colnames(datos_t) = c("semana","datos_regalos_diff")

Medias Moviles

Suavizamiento por medias moviles, para desaparecer los picos y obtener varianza constante Creamos una nueva variable suave lo añadimos como columna de datos_t de orden

datos_t$suave = ma(datos_t$datos_regalos_diff, order = 5)

Graficar datos normales y las sobre media movil dmy por el formato dia, mes, año

ggplot(data=datos_t, mapping=aes(x=dmy(semana), y=datos_regalos_diff))+
  geom_line(aes(y = datos_regalos_diff, colour = "1"))+ # linea serie original
  geom_line(aes(y = suave, colour = "2"))+ # linea medias moviles
  theme_calc()

Warning: Removed 4 row(s) containing missing values (geom_path).

Ahora tenemos una media, constante, varianza semi-constante Podemos usar un modelo arma con la siguiente serie

ggplot(data=datos_t, mapping=aes(x=dmy(semana), y=suave))+
  geom_line(aes(y = suave, colour = "2"))+ # linea medias moviles
  theme_calc()

Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Warning: Removed 4 row(s) containing missing values (geom_path).

creamor una nueva serie de tiempo y nos figamos en sus componente para determinar si es estacionaria

regalos_t_ts = ts(datos_t$suave, start=c(2016,40), frequency = 52)
plot(regalos_t_ts,xlab="Tiempo",ylab="Regalos_transformados",main="datos_transformados")

Descompanemos nuevamente esta serie Serie observada = Tendencia + Efecto est acional + Error.

regalos_t_ts_dessc = decompose(regalos_t_ts)

Graficamos

plot(regalos_t_ts_dessc, xlab='regalos_transformados')

Ajustar un modelo ARMA(p,q)

Por comodidad tomaremos los datos transformados de la diferencia solaente Buscamos los coeficientes de p y q

datos_t$datos_regalos_diff %>%
  pacf()

p = 9

datos_t$datos_regalos_diff %>%
  acf()

q = 8

Construimos el modelo ARMA(p,q)

regalos_diff_ts = datos_regalos$regalos_mexico %>% 
  diff() %>%
  ts()

regalos_diff = datos_regalos$regalos_mexico %>% 
  diff()

usamos la funcion arima para el modelo arma_p_q

m_arma_p_q = arima(regalos_diff_ts, order=c(9,0,8))

Warning in arima(regalos_diff_ts, order = c(9, 0, 8)) :
  possible convergence problem: optim gave code = 1

m_arma_p_q

Call:
arima(x = regalos_diff_ts, order = c(9, 0, 8))

Coefficients:
         ar1     ar2      ar3      ar4     ar5      ar6      ar7     ar8      ar9      ma1
      0.0245  0.1405  -0.9017  -0.3358  0.3019  -0.5348  -0.3164  0.5417  -0.2794  -0.2589
s.e.  0.0637  0.0591   0.0520   0.0650  0.0672   0.0683   0.0552  0.0557   0.0614   0.0368
          ma2     ma3     ma4      ma5     ma6     ma7      ma8  intercept
      -0.3444  0.8116  0.0009  -0.8213  0.3669  0.2355  -0.9893    -0.0364
s.e.   0.0448  0.0310  0.0503   0.0317  0.0529  0.0327   0.0404     0.0255

sigma^2 estimated as 121.5:  log likelihood = -998.7,  aic = 2035.41

Predicciones a un año con el modelo ajustado

Graficamos el modelo Predicciones de las diferencias futuras a un año

plot(forecast(m_arma_p_q,52))

Guardamos las predicciones a un año en la variable res

res <- forecast(m_arma_p_q,52)

Guardamos las esperanzas de las predicciones en un vector

forecasted = as.vector(res$mean)

base_ajustados = regalos_diff[259]+forecasted[1]# base vacia
undiff = c() # vector vacio
for (i in 1:52) {
  val = forecasted[i]
  undiff[i] = val
  base_ajustada = val
}

Datos ajustados mas predicciones

completos_ajustados=c(regalos_diff,undiff)
plot(completos_ajustados,type="l") 

Destransformaciones para visualizar las predicciones sobre la serie origina

regalos_diff = datos_regalos$regalos_mexico %>% 
  diff()

Guardamos las esperanzas de las predicciones

forecasted <- as.vector(res$mean)

las transformamos a su forma original

forecasted_original = cumsum(forecasted)-forecasted[1]+20

base_original = regalos_diff[259]+forecasted_original[1]# base vacia
undiff_original = c() # vector vacio
for (i in 1:52) {
  val_original = forecasted_original[i]
  undiff_original[i] = val_original
  base_original = val_original
}

Datos originales mas predicciones

completos_originales = c(regalos.ts,undiff_original)
plot(completos_originales,type="l")