Manejo de series de tiempo

Cargamos las librerias a usar

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

Proceso ARMA

Ejemplo 1

Simulemos el proceso $ Y_t=0.1Y_{t-1}+0.1_{t-1}+ _t $

simulcion.ma2 <- arima.sim(1000, model=list(order=c(1,0,1),ar=c(-0.1),ma=c(0.1)))
plot(simulcion.ma2)

Ahora graficamos el ACF

acf(simulcion.ma2)

Ahora, el PACF

pacf(simulcion.ma2)

Notemos que se cumple que las raíces están fuera del círculo unitario, pero no se cumple que no tengan raíces en común.

Ejemplo 2

Simulemos el proceso $ Y_t=0.9Y_{t-1}+0.4_{t-1}+ _t $

simulcion.ma2 <- arima.sim(1000, model=list(order=c(1,0,1),ar=c(-0.9),ma=c(-0.4)))
plot(simulcion.ma2)

acf(simulcion.ma2)

pacf(simulcion.ma2)

z.ar <- polyroot(c(1,0.9)) # Estacionario (Las raíces son |x|>1)
z.ma <- polyroot(c(1,-0.4)) # condicion de invertibilidad
Mod(z.ar)  
## [1] 1.111111
Mod(z.ma)
## [1] 2.5

Notemos que se cumple que las raíces están fuera del círculo unitario y no tiene raíces en común, entonces el proceso $ Y_t=0.9Y_{t-1}+0.4_{t-1}+ _t $ es estacionario.

Ejemplo 3

Datos: Serie de tiempo (1952-1988) del ingreso nacional real en China por sector (año base: 1952)

Librerias a cargar

library(AER)
## Loading required package: car
## Loading required package: carData
## Loading required package: lmtest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
data(ChinaIncome)
str(ChinaIncome)
##  Time-Series [1:37, 1:5] from 1952 to 1988: 100 102 103 112 116 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:5] "agriculture" "commerce" "construction" "industry" ...

Graficamos la serie.

transporte <- ChinaIncome[,"transport"]
ts.plot(transporte)

Parece no ser estacionario. Hacemos una transformación para tratar de confirmar la estacionariedad.

ltrans <- log(transporte)
ts.plot(ltrans)

Notamos que persiste el problema, sigue sin ser estacionario. Probemos con la diferencia:

dltrans <- diff(ltrans)
ts.plot(dltrans)
abline(h=0)

Asumamos estacionariedad (después haremos una prueba específica para verificar estacionariedad) y busquemos el mejor modelo.

Usaremos el acf y el pacf para evaluar si es MA o AR.

acf(dltrans)

pacf(dltrans)

Ajustando según la gráfica, tendríamos un proceso MA(2)

¿Qué recomienda R?

ar(dltrans)$aic
##         0         1         2         3         4         5         6         7 
##  8.140754  5.473906  1.951624  0.000000  1.990078  3.975236  5.400529  6.622877 
##         8         9        10        11        12        13        14        15 
##  6.428684  8.083145  9.975465 10.326488 12.291432 14.227677 16.183322 17.644775

Según esta recomendación, estamos ante un proceso AR(3).

modelo1 <- arima(dltrans,order = c(3,0,0))
modelo1$aic
## [1] -32.75694

Ajustemos el MA(2) y comparemos:

modelo2 <- arima(dltrans,order = c(0,0,2))
modelo2$aic
## [1] -28.01574

Recuerda: Un menor AIC es mejor. ¿Con qué modelo te quedas?

Ajustemos un MA(3).

modelo3 <- arima(dltrans,order = c(0,0,3))
print(modelo3)
## 
## Call:
## arima(x = dltrans, order = c(0, 0, 3))
## 
## Coefficients:
##          ma1      ma2      ma3  intercept
##       0.1763  -0.4596  -0.7167     0.0621
## s.e.  0.1883   0.1640   0.1925     0.0053
## 
## sigma^2 estimated as 0.01625:  log likelihood = 21.26,  aic = -32.53
modelo3$aic
## [1] -32.52547

Este modelo es mejor que el MA(2) pero peor que AR(3)

Probemos algunas combinaciones:

# Ajustando un ARMA(1,1)
modelo4 <- arima(dltrans,order = c(1,0,1))
modelo4$aic
## [1] -27.49033
# Ajustando un ARMA(2,1)
modelo5 <- arima(dltrans,order = c(2,0,1))
modelo5$aic
## [1] -32.45195
# Ajustando un ARMA(1,2)
modelo6 <- arima(dltrans,order = c(1,0,2))
modelo6$aic
## [1] -29.86264

Nos quedamos con el AR(3). Revisemos los residuos:

AR3Resid <- (modelo1$resid)
ts.plot(AR3Resid)

acf(AR3Resid)

pacf(AR3Resid)

No hay autocorrelación, no hay autocorrelación parcial.

Veamos si se trata de un ruido blanco

Box.test(AR3Resid)
## 
##  Box-Pierce test
## 
## data:  AR3Resid
## X-squared = 0.00015103, df = 1, p-value = 0.9902

¿Es ruido blanco?

Realicemos una proyección a 5 años:

pred5 <- predict(modelo1, n.ahead=5,se=T)
pred5se <- pred5$se

intervalos de confianza:

n <- (length(dltrans))
ic = pred5$pred + cbind(-2*pred5$se,2*pred5$se)
ts.plot(dltrans,pred5$pred,ic,
col=c("black","blue","red","red"))

Los intervalos son grandes, podría ser por la cantidad de datos.

A continuación se usa el modelo estimado para realizar la predicción de una sub ventana de la serie analizada:

mp <- forecast::Arima(dltrans[1:(length(dltrans)-5)], model=modelo1)
# comparamos la prediccion con lo observado:

ss <- cbind(dltrans[(length(dltrans)-4):length(dltrans)],predict(mp,n.ahead = 5)$pred)
colnames(ss) <- c("Observado","Predicción")
ss
## Time Series:
## Start = 32 
## End = 36 
## Frequency = 1 
##    Observado Predicción
## 32 0.1211453 0.08394945
## 33 0.1832397 0.05625961
## 34 0.1071942 0.05423808
## 35 0.1077345 0.06537286
## 36 0.1072015 0.07821377
ts.plot(ss,gpars=list(xlab="Año", ylab="dltrans", lty=c(1:2)))