library('ggplot2')
library('ggfortify')
library('foreign')

Ejemplo (Datos tomados del libro “Pronósticos en los negocios” de Hanke)

Cameron Consulting Corporation se especializa en ofrecer portafolios de inversión. Se desea usar la metodología de Box-Jenkins con el índice de transportación del Dow Jones. La base de datos contiene 65 promedios diarios al cierre de operaciones del mencionado índice durante los meses de verano. Los datos se encuentran el archivo Tab9-3.dta.

cameron<-read.dta("Tab9-3.dta")

Definiremos primero la serie como una serie de tiempo:

cameron<-ts(cameron[,1])

Paso 1: Identificación del Modelo

1.1. Gráfico de la Serie

El primer paso en la identificación del modelo es determinar si la serie es estacionaria, es decir, si la serie de tiempo parece variar alrededor de un nivel fijo. Es útil observar una gráfica de la serie junto con la función de autocorrelaci?n de la muestra.Examinemos el comportamiento de la serie a través del tiempo mediante un gráfico.

autoplot(cameron)

En la figura 01 aparece el gráfico de la serie original, la cual muestra un comportamiento creciente, indicando con esto la no presencia de estacionariedad. El siguiente paso en la identificación de un modelo tentativo es examinar las autocorrelaciones de muestra de datos.

1.2. Análisis de la Autocorrelación

autoplot(acf(cameron,30,main='', sub="Figura 01: Función de Autocorrelación Simple"))

El correlograma muestra un desvanecimiento lento de la función de autocorrelación indicando también no estacionariedad en la serie original.

autoplot(pacf(cameron,30,main='',sub="Figura 02: Función de Autocorrelación Parcial"))

Confirmaremos lo dicho anteriormente con la prueba de Dickey - Fuller.

1.3. Prueba de Dickey - Fuller

Para esto solicitamos el paquete fUnitRoots.

library('fUnitRoots')
## Loading required package: timeDate
## Loading required package: timeSeries
## Loading required package: fBasics

1. Planteamiento de Hipótesis

Ho: La serie es no estacionaria: Tiene raíz unitaria

H1: La serie es estacionaria: No Tiene raíz unitaria

2. Nivel de Significancia: \(\alpha\) = 0.05

3. Estadístico de Prueba

adfTest(cameron,lags=0,type=c("c"))
## Warning in adfTest(cameron, lags = 0, type = c("c")): p-value greater than
## printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: 2.2962
##   P VALUE:
##     0.99 
## 
## Description:
##  Wed Jun 24 20:48:58 2020 by user:

4. Decisión

El valor del estadístico es 2.2962, con un valor p mayor a 0.05 por lo que la hipótesis nula no se rechaza.

5. Conclusión

La serie original no es estacionaria.

Tomando Primera Diferencia

El siguiente paso es tomar la primera diferencia a la serie.

autoplot(diff(cameron))

Al tomar las primeras diferencias observamos que la serie se estabiliza en torno a un valor medio. (Figura 03)

autoplot(acf(diff(cameron)))

autoplot(pacf(diff(cameron)))

Así mismo el correlograma confirma un declinamiento rápido en el tiempo de la AFC y PAFC. (Figura 04 y 05)

Prueba de Dickey-Fuller a la serie en Primera Diferencia

1. Planteamiento de Hipótesis

Ho: La serie en primera diferencia es no estacionaria: Tiene raíz
    unitaria

H1: La serie en primera diferencia es estacionaria: No Tiene raíz   
    unitaria

2. Nivel de Significancia: \(\alpha\) = 0.05

3. Estadístico de Prueba

adfTest(diff(cameron),lags=0,type=c("c"))
## Warning in adfTest(diff(cameron), lags = 0, type = c("c")): p-value smaller than
## printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -5.8363
##   P VALUE:
##     0.01 
## 
## Description:
##  Wed Jun 24 20:49:03 2020 by user:

4. Decisión

El valor del estad?stico es -5.8363, con un valor p menor a 0.05 por lo que la hip?tesis nula se rechaza.

5. Conclusión

La serie en primera diferencia es estacionaria.

Una vez que se ha obtenido una serie estacionaria, debemos identificar la forma del modelo que se utilizará. La identificación de la forma del modelo se lleva a cabo comparando las autocorrelaciones y las autocorrelaciones parciales calculadas con los datos de la autocorrelaciones y las autocorrelaciones parciales teóricas de los diferentes modelos ARIMA.

En los correlogramas podemos observar lo siguiente: * En la ACF, la primera autocorrelación es significativa y además parece ser que los valores de la PACF disminuyen en forma oscilatoria. Esto nos lleva a pensar en un primer modelo:

\[\ y_{t}= \epsilon_{t}-\theta_{1} \epsilon_{t-1}\]

  • De igual manera, en la PACF la primera autocorrelación parcial es significativa y los valores de ACF disminuyen oscilatoriamente. Por ello proponemos un segundo modelo:

\[\ y_{t}= \phi_{0}+ \phi_{1} y_{t-1}+ \epsilon_{t}\]

Paso 2: Estimación del Modelo

Una vez que se han seleccionado los modelos tentativos, se deben estimar los par?metros para estos modelos.

Modelo 1

library('forecast')
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method                 from     
##   autoplot.Arima         ggfortify
##   autoplot.acf           ggfortify
##   autoplot.ar            ggfortify
##   autoplot.bats          ggfortify
##   autoplot.decomposed.ts ggfortify
##   autoplot.ets           ggfortify
##   autoplot.forecast      ggfortify
##   autoplot.stl           ggfortify
##   autoplot.ts            ggfortify
##   fitted.ar              ggfortify
##   fortify.ts             ggfortify
##   residuals.ar           ggfortify
library("lmtest")
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:timeSeries':
## 
##     time<-
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
arima1<-Arima(cameron, order=c(1,1,0),include.drift = TRUE)
arima1
## Series: cameron 
## ARIMA(1,1,0) with drift 
## 
## Coefficients:
##          ar1   drift
##       0.2800  1.0353
## s.e.  0.1195  0.3195
## 
## sigma^2 estimated as 3.538:  log likelihood=-130.27
## AIC=266.53   AICc=266.93   BIC=273.01
coeftest(arima1)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value Pr(>|z|)   
## ar1    0.28003    0.11946  2.3440 0.019077 * 
## drift  1.03527    0.31947  3.2406 0.001193 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo 2

arima2<-Arima(cameron, order=c(0,1,1),include.drift = TRUE)
arima2
## Series: cameron 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##          ma1   drift
##       0.2867  1.0380
## s.e.  0.1203  0.2968
## 
## sigma^2 estimated as 3.54:  log likelihood=-130.28
## AIC=266.56   AICc=266.96   BIC=273.04
coeftest(arima2)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1    0.28669    0.12029  2.3834 0.0171551 *  
## drift  1.03801    0.29678  3.4976 0.0004694 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los resultados son similares en cuanto al valor de AIC y \(\sigma^{2}\) por lo que decidimos por el modelo AR(1).

Modelo 3

arima3<-Arima(cameron,c(1,1,1),include.drift = TRUE)
arima3
## Series: cameron 
## ARIMA(1,1,1) with drift 
## 
## Coefficients:
##          ar1      ma1   drift
##       0.8356  -0.6493  1.0184
## s.e.  0.1635   0.2108  0.4659
## 
## sigma^2 estimated as 3.526:  log likelihood=-129.69
## AIC=267.38   AICc=268.06   BIC=276.01
coeftest(arima3)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1    0.83555    0.16346  5.1115 3.196e-07 ***
## ma1   -0.64926    0.21080 -3.0800   0.00207 ** 
## drift  1.01842    0.46594  2.1857   0.02884 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Paso 3: Validación del Modelo

3.1. Estacionariedad de los residuales

La siguiente figura muestra una posible estacionariedad entre los residuales del modelo. Esto lo confirmaremos con un correlograma.

checkresiduals(arima1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0) with drift
## Q* = 11.346, df = 8, p-value = 0.1828
## 
## Model df: 2.   Total lags used: 10

3.2. Autocorrelación de los residuales

autoplot(acf(arima1$residuals))

autoplot(pacf(arima1$residuals))

Los correlogramas muestran que no existe autocorrelación significativa en los residuales.

También podemos usar los paquetes ggplot2 y ggfortify para chequear la estacionariedad de los residuales.

ggtsdiag(arima(diff(cameron),c(1,0,0)))

3.2. Normalidad de los residuales

Los gráficos Q-Q son una herramienta eficaz para evaluar la normalidad. Aquí las aplicamos a los residuos.

qqnorm(arima1$residuals, sub="Figura 09: Gr?fico Q para evaluar
       normalidad");qqline(arima1$residuals)

El diagrama Q-Q de los residuos del modelo AR(1) estimado para la serie se muestra en la Figura 09. Los puntos parecen seguir la línea recta bastante de cerca, aunque parece desviars en el extremo superior. Este gráfico no nos llevar?a a rechazar la normalidad de los términos de error en este modelo.

Además, la prueba de normalidad de Shapiro-Wilk aplicada a los residuos produce un estadístico de prueba de W = 0.97148, lo que corresponde a un valor de p de 0.1445, y no rechazaría la normalidad basado en esta prueba.

shapiro.test(arima1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  arima1$residuals
## W = 0.97342, p-value = 0.1746
autoplot(arima1)

autoplot(forecast(arima1))