HANDS-ON 1
| Modelos de Supervivencia y Series de Tiempo 2022-1 |
| Ortiz Castañeda Natasha Montserrath |
ANÁLISIS DE SERIES DE TIEMPO
Extraemos una serie de tiempo de Google Trends, la cual nos dice la frecuencia con la que se buscó la palabra “pan” en México los últimos cinco años.
library(ggplot2)
library(ggthemes)
library(dplyr)
library(lubridate)
library(forecast)
library(janitor)
library(readxl)
library(dynlm)
library(stats)
setwd(file.path("C:","Users","Nath","Documents","Semestre 2022-1","series_de_tiempo","ejercicios", "data"))
pan_data<-read.csv("C://Users//Nath//Documents//Semestre 2022-1//series_de_tiempo//ejercicios//data//pandata.csv", header = TRUE)%>%
janitor::clean_names()
pan.ts <- ts(pan_data$freq_busqueda, start = c(2016,10), frequency = 48)
plot(pan.ts)
Se observa que la serie presenta tendencia aparentemente lineal y estacionalidad. Ahora analicemos su descomposición para ver si la serie puede representarse en su forma simple \(X_t = T_t + E_t + I_t +\epsilon_t\)
pan.ts.des <- decompose(pan.ts)
plot(pan.ts.des)
La tendencia puede verse creciente y decreciente entre el 2018 y pasado 2019.
shapiro.test(x=as.integer(pan.ts.des$random))##
## Shapiro-Wilk normality test
##
## data: as.integer(pan.ts.des$random)
## W = 0.92507, p-value = 6.041e-09Hecha la prueba de Shapiro, se rechaza que los errores se distribuyen normales. Por lo tanto la serie no puede modelarse como \(X_t = T_t + E_t + I_t +\epsilon_t\).
MODELO AR(p)
Para hacer un modelo autorregresivo es necesario que la serie sea estacionaria, es decir, que su media y varianza sean constantes. Para saber si nuestra serie es estacionaria, nos ayudaremos de las funciones de autocorrelación:
acf(pan.ts)
pacf(pan.ts)
El primer gráfico no luce como si se tratase de una serie estacionaria, por lo tanto tendremos que usar transformaciones que le obliguen a serlo. Nos ayudaremos de una función que nos dirá cuántos rezagos debemos hacer para que, al aplicar diferencias a la serie, tenga media constante.
ndiffs(pan.ts)## [1] 1La función nos dice que basta con un rezago:
pan.ts.dif <- diff(pan.ts, 1)
plot(pan.ts.dif)
Podemos observar que la serie ya está centrada en cero, es decir, su media ya parece ser constante. Por lo tanto, ya podemos aplicar el modelo autorregresivo.
MODELO 1: AR(1)
ar_1 <- ar(pan.ts.dif, order.max = 1)
ar_1##
## Call:
## ar(x = pan.ts.dif, order.max = 1)
##
## Coefficients:
## 1
## -0.1701
##
## Order selected 1 sigma^2 estimated as 39.15plot(forecast(ar_1, 48))
Desde aquí podemos ver que el modelo AR(1) no ajusta nada a la serie de diferencias, pero veamos cómo se ve en la serie original:
res <- forecast(ar_1, 48) #Predicción de un año
forecasted <- as.vector(res$mean)
base <- pan_data$freq_busqueda[261]
pan_data$freq_busqueda[261]+forecasted[1]## [1] 51.04307undiff <- c()
for(i in 1:48){
val <- base+forecasted[i]
undiff[i] <- val
base <- val
}
serie_completa <- c(pan_data$freq_busqueda, undiff)
plot(serie_completa, type = "l")
Hagamos uso del gráfico de autocorrelación para intentar determinar el parámetro p que ajustará mejor nuestro modelo:
pacf(pan.ts.dif)
Intentemos con \(p=3\)
MODELO 2: AR(3)
ar_3 <- ar(pan.ts.dif, order.max = 3)
ar_3##
## Call:
## ar(x = pan.ts.dif, order.max = 3)
##
## Coefficients:
## 1 2 3
## -0.2143 -0.1713 -0.1346
##
## Order selected 3 sigma^2 estimated as 37.92plot(forecast(ar_3, 48))
De igual manera se ve que no ajusta nada de nada…
MODELO 3: AR(9)
ar_9 <- ar(pan.ts.dif, order.max = 9)
ar_9##
## Call:
## ar(x = pan.ts.dif, order.max = 9)
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## -0.2549 -0.2209 -0.2045 -0.0956 -0.1069 -0.2001 -0.1492 -0.1063
## 9
## -0.1532
##
## Order selected 9 sigma^2 estimated as 36.58plot(forecast(ar_9, 48))
A pesar de que la predicción sigue siendo malísima, vemos que hay un poquito más de movimiento que en las series anteriores. Por curiosidad veamo cómo se ve en su serie original:
res <- forecast(ar_9, 48) #Predicción de un año
forecasted <- as.vector(res$mean)
base <- pan_data$freq_busqueda[261]
pan_data$freq_busqueda[261]+forecasted[1]## [1] 49.80237undiff <- c()
for(i in 1:48){
val <- base+forecasted[i]
undiff[i] <- val
base <- val
}
serie_completa <- c(pan_data$freq_busqueda, undiff)
plot(serie_completa, type = "l")
Este mal ajuste del modelo se debe a que, a pesar de haber centrado la serie, ésta no es estrictamente estacionaria como lo pide el proceso AR().
HANDS-ON 2
MODELO MA(q)
Ajustemos la serie diferenciada a un modelo MA(q).
Modelo MA(21)
Haciendo modelos de media móvil con diferentes parámetros llegamos a que el “mejor modelo” usando este proceso es con el parámetro \(q=21\).
ma_21 <- arima(pan.ts.dif, c(0,0,21))
ma_21##
## Call:
## arima(x = pan.ts.dif, order = c(0, 0, 21))
##
## Coefficients:
## ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## -0.3142 -0.1725 -0.0952 0.0030 -0.0731 -0.1882 -0.0794 -0.0004
## s.e. 0.0649 0.0667 0.0695 0.0698 0.0677 0.0729 0.0766 0.0674
## ma9 ma10 ma11 ma12 ma13 ma14 ma15 ma16
## -0.0387 0.0222 -0.0005 0.0709 0.0106 -0.0697 -0.1038 0.0523
## s.e. 0.0729 0.0764 0.0650 0.0720 0.0808 0.0893 0.0762 0.0819
## ma17 ma18 ma19 ma20 ma21 intercept
## 0.1061 0.0747 -0.1685 -0.1535 0.1181 0.0704
## s.e. 0.1042 0.0804 0.0754 0.0888 0.0776 0.0201
##
## sigma^2 estimated as 32.37: log likelihood = -823.6, aic = 1693.2plot(forecast(ma_21, 240))
Serie original:
res <- forecast(ma_21, 240) #Predicción de cinco años
forecasted <- as.vector(res$mean)
base <- pan_data$freq_busqueda[261]
pan_data$freq_busqueda[261]+forecasted[1]## [1] 51.90315undiff <- c()
for(i in 1:240){
val <- base+forecasted[i]
undiff[i] <- val
base <- val
}
serie_completa <- c(pan_data$freq_busqueda, undiff)
plot(serie_completa, type = "l")
Veamos cómo se ve un modelo ARMA(p,q):
MODELO ARMA(p,q)
ARMA(2,12)
arma_2_12 <- arima(pan.ts.dif, c(2,0,12))
arma_2_12##
## Call:
## arima(x = pan.ts.dif, order = c(2, 0, 12))
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4 ma5 ma6
## 0.5465 -0.9084 -0.8499 0.9306 -0.2962 -0.1053 -0.1645 -0.1363
## s.e. 0.0523 0.0388 0.0798 0.0932 0.1007 0.1031 0.1026 0.1042
## ma7 ma8 ma9 ma10 ma11 ma12 intercept
## -0.0047 -0.1097 -0.1219 0.0767 -0.0751 0.0807 0.0476
## s.e. 0.1001 0.0954 0.1020 0.0997 0.0866 0.0724 0.0641
##
## sigma^2 estimated as 32.78: log likelihood = -825.04, aic = 1682.09plot(forecast(arma_2_12, 240))
Modelo ARMA(9,21)
arma_9_21 <- arima(pan.ts.dif, c(9,0,21))## Warning in arima(pan.ts.dif, c(9, 0, 21)): possible convergence problem: optim
## gave code = 1arma_9_21##
## Call:
## arima(x = pan.ts.dif, order = c(9, 0, 21))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.4133 -0.0858 -0.4716 -0.6136 -0.7275 -0.6507 0.0497 -0.2601
## s.e. 0.1553 0.1336 0.1330 0.0910 0.0866 0.0733 0.1190 0.1143
## ar9 ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8
## -0.6903 0.0941 -0.1567 0.2479 0.4036 0.4984 0.0928 -0.7157 0.0217
## s.e. 0.1173 0.1719 0.1409 0.1358 0.1160 0.1207 0.1167 0.1194 0.1644
## ma9 ma10 ma11 ma12 ma13 ma14 ma15 ma16
## 0.4156 -0.5150 -0.4021 -0.2846 -0.1297 -0.1658 -0.2393 -0.0406
## s.e. 0.1620 0.1276 0.1138 0.1095 0.1026 0.1090 0.0920 0.0948
## ma17 ma18 ma19 ma20 ma21 intercept
## 0.0668 0.0107 -0.1769 -0.1436 0.1227 0.0738
## s.e. 0.0851 0.0733 0.0953 0.0832 0.1059 0.0169
##
## sigma^2 estimated as 27.73: log likelihood = -811.97, aic = 1687.94plot(forecast(arma_9_21, 240))
arma_12_21 <- arima(pan.ts.dif, c(12,0,21))## Warning in arima(pan.ts.dif, c(12, 0, 21)): possible convergence problem: optim
## gave code = 1arma_12_21##
## Call:
## arima(x = pan.ts.dif, order = c(12, 0, 21))
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): Se han producido NaNs## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.1529 0.0512 -0.3831 -0.1193 -0.2350 -0.0255 0.0063 0.0799
## s.e. NaN 0.0748 0.0663 0.0702 0.0792 0.0965 0.0854 0.0628
## ar9 ar10 ar11 ar12 ma1 ma2 ma3 ma4
## -0.4834 0.0394 -0.0577 -0.7453 -0.1861 -0.2069 0.2521 -0.0368
## s.e. 0.0643 0.0668 0.0517 0.0256 0.0581 0.0846 0.0995 0.0922
## ma5 ma6 ma7 ma8 ma9 ma10 ma11 ma12
## 0.1698 -0.3656 -0.2697 -0.0935 0.4405 -0.3221 0.0135 0.8670
## s.e. 0.0913 0.0953 0.1116 0.0564 0.0935 0.0790 0.0760 0.0631
## ma13 ma14 ma15 ma16 ma17 ma18 ma19 ma20
## -0.2575 -0.1806 -0.2901 -0.0423 0.0404 -0.2523 -0.3373 -0.1217
## s.e. 0.0836 0.0818 0.0811 0.0728 0.0751 0.0748 0.0869 0.0726
## ma21 intercept
## 0.1792 0.0734
## s.e. 0.0821 0.0170
##
## sigma^2 estimated as 26.42: log likelihood = -806.9, aic = 1683.79plot(forecast(arma_12_21, 240))
Básandonos en los valores del likelihood y aic, el último modelo es mejor que los dos primeros. En su serie original se ve de la siguiente forma:
res <- forecast(arma_12_21, 240) #Predicción de cinco años
forecasted <- as.vector(res$mean)
base <- pan_data$freq_busqueda[261]
pan_data$freq_busqueda[261]+forecasted[1]## [1] 49.73806undiff <- c()
for(i in 1:240){
val <- base+forecasted[i]
undiff[i] <- val
base <- val
}
serie_completa <- c(pan_data$freq_busqueda, undiff)
plot(serie_completa, type = "l")
De la serie original vemos que el comportamiento de los años predecidos no son muy parecidos a los años pasados. Habría que hacer más cosas para que la serie quede mejor ajustada (el hecho de que la predicción no sea buena se debe, como lo dijimos anteriormente a que la serie no es completamente estacionaria ya además tiene estacionalidad)