Cargamos las librerias a usar

library(ggplot2)
library(rmarkdown)
library(fpp2)
library(forecast)
library(timeDate)
library(gridExtra)
library(urca)

Diferenciación

Graficamos las series de tiempo correspondientes al precio de las acciones de Google y a su diferenciación

g1<-autoplot(goog200)
g2<-autoplot(diff(goog200))
grid.arrange(g1, g2, ncol = 2)

Ahora graficamos las funciones de autocorrelación

par(mfrow = c(1, 2))
acf(goog200)
acf(diff(goog200))

Se puede observar el gráfico temporal de los datos, el gráfico ACF también es útil para identificar series temporales no estacionarias. Para una serie temporal estacionaria, el ACF caerá a cero con relativa rapidez, mientras que el ACF de datos no estacionarios disminuye lentamente. Además, para datos no estacionarios, el valor de r1 suele ser grande y positiva.

Ahora procedemos a realizar las pruebas de Ljunj-Box

Box.test((goog200), lag=10, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  (goog200)
## X-squared = 1425.4, df = 10, p-value < 2.2e-16
Box.test(diff(goog200), lag=10, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  diff(goog200)
## X-squared = 11.031, df = 10, p-value = 0.3551

El ACF del precio diferenciado de las acciones de Google se parece al de una serie de ruido blanco. No hay autocorrelaciones fuera de los límites del 95% y el Ljung-Boxq∗ estadística tiene un valor p de 0,355 (para ℓ=1). Esto sugiere que el cambio diario en el precio de las acciones de Google es esencialmente una cantidad aleatoria que no está correlacionada con la de los días anteriores.

Visualización de series diferenciadas

Graficamos las ventas de antidiabetics en australia. Los registros y diferencias estacionales de los datos de ventas de A10 (antidiabético). Los logaritmos estabilizan la varianza, mientras que las diferencias estacionales eliminan la estacionalidad y la tendencia.

cbind("Sales ($million)" = a10,
      "Monthly log sales" = log(a10),
      "Annual change in log sales" = diff(log(a10),12)) %>%
  autoplot(facets=TRUE) +
  xlab("Year") + ylab("") +
  ggtitle("Antidiabetic drug sales")

Prueba de raiz unitaria

Aplicamos la prueba KPSS aprendida en clase:

goog %>% ur.kpss() %>% summary()
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 7 lags. 
## 
## Value of test-statistic is: 10.7223 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

El estadístico de prueba es mucho mayor que el valor crítico del 1%, lo que indica que se rechaza la hipótesis nula. Es decir, los datos no son estacionarios. Podemos diferenciar los datos y aplicar la prueba nuevamente.

goog %>% diff() %>% ur.kpss() %>% summary()
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 7 lags. 
## 
## Value of test-statistic is: 0.0324 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Esta vez, la estadística de prueba es pequeña y está dentro del rango que esperaríamos de datos estacionarios. Entonces podemos concluir que los datos diferenciados son estacionarios.

Este proceso de utilizar una secuencia de pruebas KPSS para determinar el número apropiado de primeras diferencias lo lleva a cabo la función ndiffs().

ndiffs(goog)
## [1] 1

Modelo ARIMA

Graficamos los cambios porcentuales trimestrales en el gasto de consumo de Estados Unidos. Aunque se trata de una serie trimestral, no parece haber un patrón estacional, por lo que ajustaremos un modelo ARIMA no estacional.

autoplot(uschange[,"Consumption"]) +
  xlab("Year") + ylab("Quarterly percentage change")

Vemos si sugiere diferenciacion

ndiffs(uschange)
## [1] 0

No sugiere diferenciacion

Seleccionamos un modelo automaticamente con la función “autoarima”

fit <- auto.arima(uschange[,"Consumption"], seasonal=FALSE)
fit
## Series: uschange[, "Consumption"] 
## ARIMA(1,0,3) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     ma2     ma3    mean
##       0.5885  -0.3528  0.0846  0.1739  0.7454
## s.e.  0.1541   0.1658  0.0818  0.0843  0.0930
## 
## sigma^2 = 0.3499:  log likelihood = -164.81
## AIC=341.61   AICc=342.08   BIC=361

Nos indica que es un modelo ARIMA (1,0,3)

\[ y_t=C+0.589y_{t-1}-0.353\varepsilon_{t-1}+0.0846\varepsilon_{t-2}+0,174\varepsilon_{t-3}+\varepsilon_t \] donde \(C=0.745*(1-0.589)=0.307\) y \(\varepsilon_t\) es ruido blanco con una desviación estándar de \(0,592=\sqrt{0.307}\)

Luego, mostramos el pronostico del mismo

fit %>% forecast(h=10) %>% autoplot(include=80)

Pero, es este el mejor modelo? vamos a analizarlo, para esto, seleccionaremos el mejor modelo manualmente.

Graficamos las funciones de autocorrelacion y autocorrelacion parcial

par(mfrow = c(1, 2))
acf(uschange[,"Consumption"])
pacf(uschange[,"Consumption"])

Vemos que hay tres picos en el ACF, seguidos por un pico casi significativo en el retraso 4. En el PACF, hay tres picos significativos, y luego no hay picos significativos (aparte de uno justo fuera de los límites en retraso 22). Podemos ignorar un pico significativo en cada gráfico si está justo fuera de los límites y no en los primeros rezagos. Después de todo, la probabilidad de que un pico sea significativo por casualidad es de aproximadamente uno entre veinte, y estamos trazando 22 picos en cada gráfico. El patrón en los primeros tres picos es el que esperaríamos de ARIMA(3,0,0), ya que el PACF tiende a disminuir. Entonces, en este caso, ACF y PACF nos llevan a pensar que un modelo ARIMA(3,0,0) podría ser apropiado.

(fit2 <- Arima(uschange[,"Consumption"], order=c(3,0,0)))
## Series: uschange[, "Consumption"] 
## ARIMA(3,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2     ar3    mean
##       0.2274  0.1604  0.2027  0.7449
## s.e.  0.0713  0.0723  0.0712  0.1029
## 
## sigma^2 = 0.3494:  log likelihood = -165.17
## AIC=340.34   AICc=340.67   BIC=356.5

En realidad, este modelo es ligeramente mejor que el modelo identificado por auto.arima()(con un valor AICc de 340,67 en comparación con 342,08). La función auto.arima() no encontró este modelo porque no considera todos los modelos posibles en su búsqueda. Puedes hacer que funcione más usando los argumentos stepwise=FALSEy approximation=FALSE:

(fit3 <- auto.arima(uschange[,"Consumption"], seasonal=FALSE,
                    stepwise=FALSE, approximation=FALSE))
## Series: uschange[, "Consumption"] 
## ARIMA(3,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2     ar3    mean
##       0.2274  0.1604  0.2027  0.7449
## s.e.  0.0713  0.0723  0.0712  0.1029
## 
## sigma^2 = 0.3494:  log likelihood = -165.17
## AIC=340.34   AICc=340.67   BIC=356.5

Esta vez, con auto.arima() encontramos el mismo modelo que supusimos a partir de los gráficos ACF y PACF

forecast1 <- fit %>% forecast(h=10) %>% autoplot(include=50)
forecast2 <- fit3 %>% forecast(h=10) %>% autoplot(include=50)
grid.arrange(forecast1, forecast2, ncol = 1)