Taller 6 - Series de tiempo - Modelo Arima

Carga librerías

setwd("E:/Personal/especializacion/series_de_tiempo/taller6/")

#library(fma)
library(fpp2)
library(gridExtra)
library(grid)
library(ggplot2)
library(lattice)
library(tseries)
library(urca)

PUNTO 3

Para la siguiente serie, encuentre una transformación de Box-Cox y un orden de diferenciación apropiados para obtener datos estacionarios.

a.usneletelec

b.enplanements

c.visitors

ANALISIS SERIE USNELETELEC

Descripción de serie

#?usnetelec
#usnetelec
ggtsdisplay(usnetelec)

La serie describe generación de electricidad anual en kwh PARA LOS AÑOS 1949 A 2003

La serie de tiempo no presenta estacionalidad, no es periódica, no es necesaria transformación box-cox, puesto que se observa una varianza homogénea

se observa una tendencia creciente casi lineal (la serie no es estacionaria en media) se requiere diferenciar , por tanto se puede decir que la serie no es estacionaria.

El diagrama ACF confirma que la serie no es estacionaria

Utilizamos las funciones ndiffs, que calculan cada una el número de diferenciaciones regulares que se necesitan llevar a cabo para que la serie sea estacionaria.

Identificación de numero de diferenciaciones regulares

ndiffs(usnetelec)

## [1] 1

El calculo indica que la serie necesita una(1) diferenciación regular, no es necesaria una diferenciación estacional puesto que la serie no lo es.

Diferenciación

#diferenciando
diff.co2ts<-autoplot(diff(usnetelec) )
diff.co2ts

autoplot(acf(diff(usnetelec), plot = FALSE))

Se observa eliminación en el componente de tendencia

Prueba de estacionariedad

Una vez eliminado la componente de tendencia , observamos que esta serie se parece bastante a una serie estacionaria, ya que parece ser constante en media y varianza, para asegurarnos aplicamos un test de estacionariedad, el test de KPSS (Kwiatkowski-Phillips-Schmidt-Shin).

Para confirmar usamos test kpss

diff(usnetelec) %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.1585 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

\(H_o\): La serie es estacionaria

\(H_a\): La serie es no estacionaria

Como el valor estadístico calculado |0.1585| en valor absoluto es menor que los valores críticos, acepto la hipótesis nula. Por lo tanto la serie es estacionaria con una diferenciación.

ANALISIS SERIE ENPLANEMENTS

Descripción de la serie

#enplanements
ggtsdisplay(enplanements)

Se observa una tendencia ascendente (la serie no es estacionaria en media), por tanto no estacionaria

Se observa una varianza no homogénea, las variaciones son grandes a medida que los años van pasando, aplicamos por tanto transformación box-cox

También podemos apreciar que existe la componente estacional

Componentes de la serie enplanements

Se observa tendencia, estacionalidad y varianza no homogénea

plot(decompose(enplanements))

Identificación de diferenciaciones a realizar

Tenemos que eliminar la tendencia y la estacionalidad de la serie para que pueda ser estacionaria, diferenciando la serie lograremos que sea convertida a estacionaria, usamos las funciones ndiffs y nsdiffs para calcular diferenciaciones en tendencia y estacionales respectivamente

ndiffs(enplanements)

## [1] 1

nsdiffs(enplanements)

## [1] 1

frequency(enplanements)

## [1] 12

Los resultados dicen que usaremos una diferenciación estacional y una primera diferencia, Aplicamos box-cox a la serie para estabilizar la varianza

Aplicación Box-Cox, diferencia estacional y diferenciación regular

enpl.lambda <- BoxCox.lambda(enplanements)
enpl.lambda

## [1] -0.2269461

enpl.estacionaria <- enplanements %>% BoxCox(enpl.lambda) %>% diff(lag=12) %>% diff(1)
ggtsdisplay(enpl.estacionaria)

El lambda de la transformación cox-box es \(\lambda=-0.227\)

Se ha eliminado la tendencia y se estabiliza la varianza , para confirmar la estacionariedad usaremos el test kpss

Test KPSS para confirmar estacionariedad

enpl.estacionaria %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0424 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

El test estadístico es |0,0424| es menor que los valores críticos, acepto la hipótesis nula. Por lo tanto la serie ya es estacionaria con una diferenciación.

ANALISIS SERIE VISITORS

Descripción de la serie

La serie visitors es la serie conformada por visitantes mensuales australianos a corto plazo en el extranjero. Mayo de 1985-abril de 2005

#?visitors
ggtsdisplay(visitors)

Se observa una tendencia ascendente (la serie no es estacionaria en media), por tanto no estacionaria, observamos varianza no homogénea, las variaciones van aumentando a medida que los años van pasando, debemos aplicar transformación para estabilizar la varianza y lo haremos con box-cox. Se observa unos picos cada 12 lags.

Identificación de componentes de la serie

Se observa los componentes tendencia, estacionalidad en la serie

plot(decompose(visitors))

Identificación de diferenciaciones regulares y de estacionalidad a realizar

Tenemos que eliminar la tendencia y la estacionalidad de la serie para que pueda ser estacionaria, diferenciando la serie lograremos que sea convertida a estacionaria, usamos las funciones ndiffs y nsdiffs para calcular diferenciaciones en tendencia y estacionales respectivamente

ndiffs(visitors)

## [1] 1

nsdiffs(visitors)

## [1] 1

Los resultados dicen que usaremos una diferenciación regular y otra estacional

Aplicación de transformación y diferenciación

visitors.lambda <- BoxCox.lambda(visitors)
visitors.estacionaria <- visitors %>% BoxCox(visitors.lambda) %>% diff(12) %>% diff(1)
ggtsdisplay(visitors.estacionaria )

Se ha eliminado la tendencia y se estabiliza la varianza , para confirmar la estacionariedad usaremos el test kpss

Test de estacionariedad

visitors.estacionaria %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 0.0158 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

El test estadístico es |0,0158| es menor que los valores críticos, acepto la hipótesis nula. Por lo tanto la serie ya es estacionaria con una diferenciación regular y una diferenciación estacional.

PUNTO 4

Para los datos de enplanements, escriba las diferencias que eligió anteriormente utilizando la notación de operador de retroceso.

La notación de retroceso es común para representar modelos ARIMA, el cual usa el operador B el cual desplaza datos hacia atrás un periodo.

Una diferencia de orden d puede ser escrita como

\((1-B)(1-B^m)y_{t}=(1-B-B^m+B^{m+1})y_{t}\)

La serie enplanements necesitó

• una primera diferencia regular (d=1)

• una diferencia estacional (m=12)

• una transformación box-cox

Para la serie enplanements tenemos

\((1-B)(1-B^{12})y_{t}=(1-B-B^{12}+B^{13})y_{t}=e_{t}\)

Aunque creemos que la ecuación correcta es:

\((1-B)(1-B^{12})( \frac{y_{t}^{-0.227}-1}{-0.227} )=(1-B-B^{12}+B^{13})( \frac{y_{t}^{-0.227}-1}{-0.227} )=e_{t}\)

Puesto que tiene una transformación cox-box con \(\lambda=-0.227\) Profesor esto nos gustaría nos revise si esta bien

PUNTO 5

Para sus datos minoristas, encuentre el orden apropiado de diferenciación (después de la transformación si es necesario) para obtener datos estacionarios.

Lectura y descripción de la serie

Se realiza la lectura de los datos de la base de datos base2.xlsx y escogemos los datos de la columna A3349399C referente a venta de ropa al por menor

library(fpp2)
retaildata<-readxl::read_excel("base2.xlsx",skip=1)
#retaildata
miSerieRopa <- ts(retaildata[,"A3349399C"], frequency=12, start=c(1982,4))
plot(decompose(miSerieRopa))

ggtsdisplay(miSerieRopa)

Se observa tendencia ascendente (la serie no es estacionaria en la media), por tanto no estacionaria , los datos presentan estacionalidad (retardo significativo cada 12 lags ) Se observa que la varianza no es homogénea,la varianza incrementa conforme pasan los años, por tanto se aplicara transformación box-cox, se eliminará la tendencia y la estacionalidad de la serie para que pueda ser estacionaria, diferenciando la serie, usamos las funciones ndiffs y nsdiffs para calcular diferenciaciones en tendencia y estacionales respectivamente

Identificación de cantidad de diferenciaciones regulares y estacionales

ndiffs(miSerieRopa)

## [1] 1

nsdiffs(miSerieRopa)

## [1] 1

Los resultados dicen que usaremos una primera diferenciación y otra estacional

Por tanto la serie estacionaria se obtiene de la siguiente forma

miSerieRopa.lambda <- BoxCox.lambda(miSerieRopa)
miSerieRopa.lambda

## [1] 0.02074707

miSerieRopa.estacionaria <- miSerieRopa %>% BoxCox(miSerieRopa.lambda) %>% diff(lag=12) %>% diff(1)
ggtsdisplay(miSerieRopa.estacionaria )

El mejor lambda box-cox para esta serie es 0.02074707

Se observa que se elimina la tendencia y se estabiliza la varianza , para confirmar la estacionariedad usaremos el test kpss

Test de estacionariedad

miSerieRopa.estacionaria %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0316 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

El test estadístico es |0.0316 | es menor que los valores críticos, acepto la hipótesis nula. Por lo tanto la serie ya es estacionaria con una transformación , una diferenciación regular y una diferenciación estacional

PUNTO 6

Considere wmurders la cantidad de mujeres asesinadas cada año (por cada 100,000 habitantes estándar) en los Estados Unidos.

a. Al estudiar los gráficos apropiados de la serie en R, encuentre un ARIMA apropiado.

b. ¿Deberías incluir una constante en el modelo? Explique.

c. Escriba este modelo en términos del operador de rezagos.

d. Ajuste el modelo usando R y examine los residuos. ¿El modelo es satisfactorio?

e. Pronostica las siguientes tres observaciones.

f. Grafique los pronósticos

g. ¿auto.arima() Da el mismo modelo que has elegido? Si no, ¿qué modelo crees que es mejor?

Descripción de la serie

Se observa una tendencia creciente hasta 1975, luego se estabiliza y después decrece en la década 1990 hasta 2000, se observa por tanto que no es una serie estacionaria.

No hay incremento en varianza a través de los años, por tanto no es necesaria una transformación box-cox, no se observa estacionalidad.

a. Al estudiar los gráficos apropiados de la serie en R, encuentre un ARIMA apropiado.

#plot(decompose(wmurders))
ggtsdisplay(wmurders)

Se eliminará la tendencia de la serie para que pueda ser estacionaria, diferenciando la serie, usamos las funciones ndiffs para calcular diferenciaciones en tendencia

Identificación de diferenciaciones a realizar

ndiffs(wmurders)

## [1] 2

Los resultados recomiendan que diferencie dos veces (d=2) para que sea estacionaria

Aplicación de diferenciación

Por tanto la serie estacionaria se obtiene de la siguiente forma

wmurders.estacionaria <- wmurders %>% diff(1) %>% diff(1)
ggtsdisplay(wmurders.estacionaria)

Se observa que se elimina la tendencia, para confirmar la estacionariedad usaremos el test kpss

Test de estacionaridad

wmurders.estacionaria %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.0458 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

El test estadístico es |0.0458 | es muy pequeño, menor que los valores críticos, acepto la hipótesis nula. Por lo tanto la serie ya es estacionaria con dos diferenciaciones.

Obtención de Arima apropiado serie wmurders, ajuste en r y observación de residuos

Se observa después de la diferenciación en el ACF un comportamiento sinusoidal por tanto el orden lo da el PACF, por lo cual nos inclinamos por un AR con p=1 (tenemos 1 lag significativo) y como tenemos d=2 (dos diferenciaciones) un primer modelo puede ser

ARIMA(1,2,0)

arima.murders.1 <- Arima( wmurders, order=c(1,2,0) )
arima.murders.1

## Series: wmurders 
## ARIMA(1,2,0) 
## 
## Coefficients:
##           ar1
##       -0.6719
## s.e.   0.0981
## 
## sigma^2 estimated as 0.05471:  log likelihood=2
## AIC=0   AICc=0.24   BIC=3.94

Observamos también un pico significativo en el ACF ( q=1 ) en el MA , por tanto un segundo modelo puede ser:

ARIMA(1,2,1)

arima.murders.2 <- Arima( wmurders, order=c(1,2,1) )
arima.murders.2

## Series: wmurders 
## ARIMA(1,2,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.2434  -0.8261
## s.e.   0.1553   0.1143
## 
## sigma^2 estimated as 0.04632:  log likelihood=6.44
## AIC=-6.88   AICc=-6.39   BIC=-0.97

Observación de Residuos de los dos modelos a estudiar

Residuos obtenidos para modelo ARIMA(1,2,0)

checkresiduals( arima.murders.1 )

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,2,0)
## Q* = 16.452, df = 9, p-value = 0.05802
## 
## Model df: 1.   Total lags used: 10

Residuos obtenidos para modelo ARIMA(1,2,1)

checkresiduals( arima.murders.2 )

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

Los residuales presentan comportamiento ruido blanco en los dos modelos, la prueba de Ljong box , p-value > 0.05 se acepta la hipótesis nula de que los datos se distribuyen de manera independiente, no hay correlación en los residuos, por tanto los datos son ruido blanco, En general los residuos de los dos modelos presentan comportamiento deseado y presenta normalidad

Ahora bien, comparando AIC de los dos modelos obtenemos \[AIC\_arima.murders.2 < AIC\_arima.murders.1\]

Por tanto el mejor modelo encontrado en estos dos casos es el generado por ARIMA(1,2,1)

Inclusión de una constante en el modelo

No se debe incluir constante, investigando encontramos que Si d = 0, entonces se incluye la constante c; si d>=1 entonces la constante c se pone a cero. Las variaciones en el modelo actual se consideran variando p y / o q del modelo actual en $ 1$ e incluyendo / excluyendo c del modelo actual.

Escriba este modelo en términos del operador de rezagos.

Usaremos el modelo Arima(1,2,1), la ecuación queda como en la figura

\((1-\phi_{1}B)(1-B)^2y_{t}=c +(1+\theta_1B)\varepsilon_{t}\)

Puesto que tenemos un AR de orden 1, dos diferenciaciones y un MA de orden 1

Pronóstico de las siguientes tres observaciones.

Para nuestro modelo Arima(1,2,1) escogido obtenemos las siguientes tres observaciones

forecast(arima.murders.2 , h=3)

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005       2.470660 2.194836 2.746484 2.048824 2.892496
## 2006       2.363106 1.986351 2.739862 1.786908 2.939304
## 2007       2.252833 1.765391 2.740276 1.507354 2.998313

Grafico los pronósticos

autoplot(forecast(arima.murders.2 , h=3))

Observamos en la grafica una tendencia a que la cantidad de mujeres asesinadas cada año (por cada 100,000 habitantes estándar) en los Estados Unidos esta bajando en la ultima decada y lo confirman los pronosticos de los siguientes tres años.

¿auto.arima() Da el mismo modelo que has elegido? Si no, ¿qué modelo crees que es mejor?**

Este tercer modelo lo realizaremos con auto.arima de R obtenemos

arima.murders.3 <- auto.arima( wmurders, stepwise = FALSE, approximation = FALSE)
arima.murders.3

## Series: wmurders 
## ARIMA(0,2,3) 
## 
## Coefficients:
##           ma1     ma2      ma3
##       -1.0154  0.4324  -0.3217
## s.e.   0.1282  0.2278   0.1737
## 
## sigma^2 estimated as 0.04475:  log likelihood=7.77
## AIC=-7.54   AICc=-6.7   BIC=0.35

El modelo encontrado con autoarima con stepwise= FALSO recomienda un ARIMA (0,2,3) obtenemos un AIC=-7.54, dando un resultado mejor que los dos modelos que estamos trabajando ARIMA(1,2,1) y ARIMA(1,2,0)

Los residuales del autoarima , se muestran en pantalla

checkresiduals( arima.murders.3 )

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,2,3)
## Q* = 10.706, df = 7, p-value = 0.152
## 
## Model df: 3.   Total lags used: 10

Se observa unos residuales presentando ruido blanco, la prueba de Ljong box , p-value > 0.05 se acepta la hipótesis nula de que los datos se distribuyen de manera independiente, no hay correlación en los residuos, por tanto los datos son ruido blanco, En general los residuos de los dos modelos presentan comportamiento deseado y presenta normalidad

PUNTO 8

La producción anual de carbón bituminoso en los Estados Unidos entre 1920 y 1968 se encuentra en el conjunto de datos bicoal.

a. Producir un diagrama de tiempo de los datos.

b. Decide ajustar el siguiente modelo a la serie:

\(y_{t}=C+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\phi_{3}y_{t-3}+\phi_{4}y_{t-4}+\varepsilon_{t}\)

c. Explique por qué se eligió este modelo utilizando ACF y PACF

d. Los últimos cinco valores de la serie se dan a continuación.

AÑO	1964	1965	1966	1967	1968
Millones de Toneladas	467	512	534	552	545

Los parámetros estimados son

\(C = 162, \phi_{1} = 0.83, \phi_{2} = -0.34, \phi_{3} = 0.55, \phi_{4} = -0.38.\)

Sin usar la función forecast, calcule pronósticos para los próximos tres años (1969-1971).

e. Ahora ajuste el modelo en R y obtenga los pronósticos del mismo modelo. ¿Cómo son diferentes de los tuyos?

Diagrama de Tiempos y descripción

  ggtsdisplay(bicoal)

No se observa un patrón en la serie, podría ser estacionaria, para efectos del punto vamos a aplicar un test KPSS de estacionariedad para confirmar esto.

bicoal %>% ur.kpss() %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.0905 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Efectivamente elvalor del test estadístico (0.0905) es muy bajo comparado con los valores críticos, esto significa que la serie es estacionaria.

Modelo Arima relacionado con \(y_t\)

\(\varepsilon_{t}\) es ruido blanco, se observa un modelo autoregresivo de orden p=4

AR(4) = ARIMA(4,0,0)

Explique por qué se eligió este modelo utilizando ACF y PACF

El ACF decae exponencialmente (podriamos decir que el ACF en conjunto es sinusoidal), el orden lo da EL PACF , el último pico significativo en el PACF es el 4 y no hay ninguno mas alla de este pico , nos inclinamos por un AR, obteniendo el modelo del punto b ARIMA(4,0,0)

Pronóstico de los próximos cuatro años manualmente

Para hacer los cálculos creamos una función con la ecuación de pronóstico

\(y_{t}=C+\phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\phi_{3}y_{t-3}+\phi_{4}y_{t-4}+\varepsilon_{t}\)

yt.bicoal <- function ( yt_1, yt_2, yt_3, yt_4){
  C=162
  phi1 = 0.83
  phi2 = -0.34
  phi3 = 0.55
  phi4 = -0.38
  y <- C+ phi1*yt_1 + phi2*yt_2 + phi3 *yt_3 + phi4*yt_4
  return (y)
}

y1965 = 512
y1966 = 534
y1967 = 552
y1968 = 545

y1969 = yt.bicoal(y1968, y1967, y1966, y1965);y1969

## [1] 525.81

y1970 = yt.bicoal(y1969, y1968, y1967, y1966);y1970

## [1] 513.8023

y1971 = yt.bicoal(y1970, y1969, y1968, y1967);y1971

## [1] 499.6705

y1972 = yt.bicoal(y1971, y1970, y1969, y1968);y1972

## [1] 484.1292

Los resultados para los siguientes cuatro años con el modelo del punto d se resumen en la siguiente tabla

AÑO	1969	1970	1971	1972
Millones de Toneladas	525.81	513.80	499.67	484.13

Ajuste del modelo en R :

arima.bicoal <- Arima( bicoal, order=c(4,0,0) )
arima.bicoal

## Series: bicoal 
## ARIMA(4,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      mean
##       0.8334  -0.3443  0.5525  -0.3780  481.5221
## s.e.  0.1366   0.1752  0.1733   0.1414   21.0591
## 
## sigma^2 estimated as 2795:  log likelihood=-262.05
## AIC=536.1   AICc=538.1   BIC=547.45

Se observa los mismos valores para \(\phi_{1}, \phi_{2}, \phi_{3}, \phi_{4}\), con \(\mu = 481.5221\)

valor de constante

\(C =\mu * (1-\phi_{1}-\phi_{2}-\phi_{3}-\phi_{4})\)

\(C =481.5221 * (1-0.8334+0.3443-0.5525+0.3780)= 161.984\)

Realizamos pronóstico para los próximos 4 años usando forecast

autoplot(forecast(arima.murders.2 , h=4))

pronost.bicoal<-forecast(arima.bicoal, h=4)
pronost.bicoal

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1969       527.6291 459.8804 595.3779 424.0164 631.2419
## 1970       517.1923 429.0014 605.3832 382.3160 652.0686
## 1971       503.8051 412.4786 595.1315 364.1334 643.4768
## 1972       489.2909 390.4658 588.1160 338.1510 640.4309

Los pronósticos son muy semejantes, de pronto puede ser tema de redondeos para valores de \(\phi\) y constante C en general se observa una tendencia a la baja en la producción anual de carbón bituminoso

!Muchas gracias Profesor¡