Basado en las notas de Luis Rincón.

Paquetería requerida • stats: para el análisis de series de tiempo • tseries: para pruebas de hipótesis acerca de la estacionariedad • nortest: para la prueba de hipótesis de normalidad • FinTS: para la prueba de hipótesis para homocedasticidad

library("stats")
library("tseries")
## Warning: package 'tseries' was built under R version 4.4.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library("nortest")
library("FinTS")
## Warning: package 'FinTS' was built under R version 4.4.3
## Cargando paquete requerido: zoo
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library('forecast')
## Warning: package 'forecast' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'forecast'
## The following object is masked from 'package:FinTS':
## 
##     Acf

Usaremos datos referentes a el número total de pasajeros mensuales de aerolíneas internacionales de 1949 a 1960:

data("AirPassengers")
View(AirPassengers)

Los datos ya están en formato de Serie de Tiempo:

serie <- AirPassengers
class(serie)
## [1] "ts"

La función ts( ), convierte a formato serie de tiempo.

serie
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
plot(serie, col = "green", lty = 1, lwd = 4, ylab = "Pasajeros",
xlab = "Tiempo en meses", legend = NULL, main="Serie de Tiempo Inicial/Original")

Eliminación de Tendencia y Estacionalidad.

La función ‘stl( )’, puede hacer una descomposición que considera estacionalidad:

#stl(serie,s.window="periodic")
stl(serie,s.window="periodic")$time.series[1:10,]
##         seasonal    trend   remainder
##  [1,] -25.497718 127.1873  10.3103699
##  [2,] -35.220935 126.6495  26.5714028
##  [3,]  -3.027478 126.1117   8.9157625
##  [4,]  -8.299054 126.1989  11.1001224
##  [5,]  -5.737289 126.2861   0.4511408
##  [6,]  32.336634 126.7330 -24.0696572
##  [7,]  70.243882 127.1799 -49.4237797
##  [8,]  68.049433 127.4162 -47.4656492
##  [9,]  17.438326 127.6525  -9.0908608
## [10,] -21.063432 129.0186  11.0448425
?stl
## starting httpd help server ... done
plot(stl(serie,s.window="periodic"),labels=c("Datos","Estacion.","Tenden.","Residuo"))

Que parece ser un poco más estacionaria.

Podemos hacer diferencias para eliminar tendencias y estacionalidad. La función diff( ) sirve para obtener la serie diferenciada, lag indica el retraso (por ejemplo para un modelo con estacionalidad s podríamos usar lag=s) y differences indica el orden de las diferencias, es decir, el número de veces que se está diferenciando. Obtendremos la serie con lag de 12 para eliminar estacionalidad:

S=12, En un modelo SARIMA

diflag12<-diff(serie, lag = 12, differences = 1)
plot(diflag12, col = "darkred", lty = 1, lwd = 4, ylab = "Diferencias",
xlab = "Tiempo en meses", legend = NULL, main="Serie, ¿sin Estacionalidad?")
grid(10, 12, lwd = 2)

diflag1<-diff(diflag12, lag = 1, differences = 1)
plot(diflag1, col = "darkred", lty = 1, lwd = 4, ylab = "Diferencias",
xlab = "Tiempo en meses", legend = NULL, main="Serie, ¿sin Tendencia?")
grid(10, 12, lwd = 2)

Si se hacen diferencias 2 veces sobre la serie sin estacionalidad “difflag12” entonces se usa la instrucción siguiente:

diflag2<-diff(diflag12, lag = 1, differences = 2)
plot(diflag2, col = "darkred", lty = 1, lwd = 4, ylab = "Diferencias",
xlab = "Tiempo en meses", legend = NULL, main="Serie, sin Tendencia")

Las decisiones para eliminar la tendencia y estacionalidad para obtener una serie estacionaria de manera gráfica son subjetivas, sin embargo, existen pruebas de hipótesis estadísticas que nos ayudan a tomar dichas decisiones. Con las pruebas de raíz unitaria Dickey-Fuller “DF”, Augmented Dickey-Fuller “ADF” y PhillipsPerron “PP”, podemos probar si la serie diflag1 es estacionaria (H0: la serie no es estacionaria vs H1: la serie es estacionaria):

\[ H_0: Hay raiz unitaria \] \[ H_a: No hay raiz unitaria \]

# Prueba DF:
adf.test(diflag1, k=0)
## Warning in adf.test(diflag1, k = 0): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diflag1
## Dickey-Fuller = -15.558, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary
# Se rechaza que la serie no es estacionaria
# Prueba DF aumentada:

Tendencia d=1, S=12.

Selección del Modelo y Estimación de Parámetros.

Graficaremos las autocorrelaciones y autocorrelaciones parciales para seleccionar los valores de los parámetros en nuestro modelo ARMA (p,q):

serieamodelar<-scale(diflag1, scale = FALSE)
serieamodelar[1:15]
##  [1]  4.8167939  0.8167939 -3.1832061 -2.1832061  9.8167939  7.8167939
##  [7] -0.1832061 -0.1832061 -8.1832061 -4.1832061 11.8167939  7.8167939
## [13] -6.1832061 12.8167939 -9.1832061
plot(serieamodelar, col = "darkred", lty = 1, lwd = 4, ylab = "Serie Diferenciada",
xlab = "Tiempo en meses", legend = NULL, main="Serie Estacionaria", type = "l")

# MA(q)
# Parámetro q
acf(ts(serieamodelar,freq=1), lag.max=30, main="ACF del proceso")

MA(1)

# Parámetro p
pacf(ts(serieamodelar,freq=1), lag.max=30, main="PACF del proceso")

AR(1)

De las autocorrelaciones se sale la ACF para el lag 1 y 23, se sugiere al menos MA(1) y en cuanto a las autocorrelaciones parciales PACF se sugiere también AR de al menos 1, pues los lags 20 y 22 también salen de las bandas, haremos un ajuste para un ARMA(1,1) por máxima verosimilitud:

En realidad SARIMA12. Notemos que la serie a modelar ya no tiene tendencia pues se diferencio, por ese motivo el modeo a ajustar es el modelo ARMA(1,1)

ARMA11<-arima(serieamodelar,order=c(1,0,1), include.mean= FALSE)
ARMA11
## 
## Call:
## arima(x = serieamodelar, order = c(1, 0, 1), include.mean = FALSE)
## 
## Coefficients:
##           ar1      ma1
##       -0.2957  -0.0133
## s.e.   0.3935   0.4238
## 
## sigma^2 estimated as 137:  log likelihood = -508.17,  aic = 1022.34

\[H_0: \phi_1 = 0\] \[H_a: \phi_1 \neq 0\]

-0.2957/0.3935
## [1] -0.7514612
2*pnorm(0.3935,lower.tail = FALSE)
## [1] 0.6939502

No se rechaza H0

\[H_0: \theta_1 = 0\] \[H_a: \theta_1 \neq 0\]

-0.0133/0.4238
## [1] -0.03138273
2*pnorm(-0.03138273,lower.tail = T)
## [1] 0.9749643

Como en la serie original lo que hicimos fue diferenciar una vez de acuerdo a la estacionalidad, o sea con lag de 12, y luego una vez para eliminar tendencia usando un lag de 1, en realidad pudimos haber introducido el modelo con los datos orginales y obtener el mismo modelo ARMA(1,1) visto introduciendo al modelo como un SARIMA((1,1,1),(0,1,0)), donde el primer paréntesis indica que es un ARMA diferenciando una vez (el segundo 1 del paréntesis indica esto). El segundo paréntesis indica que en la parte estacional se diferenció una vez (esta es la diferencia con lag 12).

ARMA11forma2<-arima(serie, order=c(1,1,1),
seasonal = list(order = c(0, 1, 0), period = NA),
include.mean= FALSE)
ARMA11forma2
## 
## Call:
## arima(x = serie, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), period = NA), 
##     include.mean = FALSE)
## 
## Coefficients:
##           ar1      ma1
##       -0.3009  -0.0073
## s.e.   0.3835   0.4133
## 
## sigma^2 estimated as 137:  log likelihood = -508.2,  aic = 1022.39
# Notemos que da lo mismo si incluyo o no la media arriba
ARMA11forma2<-arima(serie, order=c(1,1,1),
seasonal = list(order = c(0, 1, 0), period = NA),
include.mean= TRUE)
ARMA11forma2
## 
## Call:
## arima(x = serie, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), period = NA), 
##     include.mean = TRUE)
## 
## Coefficients:
##           ar1      ma1
##       -0.3009  -0.0073
## s.e.   0.3835   0.4133
## 
## sigma^2 estimated as 137:  log likelihood = -508.2,  aic = 1022.39

No es aconsejable hacer uso del AUTOARIMA unicamente. Se aconseja que se haga uso de varios métodos para posteriormente comparar

mod2<-auto.arima(serie)
mod2
## Series: serie 
## ARIMA(2,1,1)(0,1,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.5960  0.2143  -0.9819
## s.e.  0.0888  0.0880   0.0292
## 
## sigma^2 = 132.3:  log likelihood = -504.92
## AIC=1017.85   AICc=1018.17   BIC=1029.35
ARMA11.pred<-predict(ARMA11forma2,n.ahead=12)
ARMA11.pred
## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1961 444.3670 418.2566 446.2898 488.2798 499.2828 562.2819 649.2822 633.2821
##           Sep      Oct      Nov      Dec
## 1961 535.2821 488.2821 417.2821 459.2821
## 
## $se
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1961 11.70537 14.23287 16.93806 19.11351 21.10701 22.91630 24.59607 26.16732
##           Sep      Oct      Nov      Dec
## 1961 27.64968 29.05644 30.39819 31.68316
AUARMA11.pred<-predict(mod2,n.ahead=12)
AUARMA11.pred
## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1961 445.6349 420.3950 449.1983 491.8399 503.3945 566.8624 654.2602 638.5975
##           Sep      Oct      Nov      Dec
## 1961 540.8837 494.1266 423.3327 465.5076
## 
## $se
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1961 11.50524 13.50261 15.15798 16.24046 17.04072 17.63388 18.08602 18.43632
##           Sep      Oct      Nov      Dec
## 1961 18.71242 18.93348 19.11322 19.26159

Verificación de Supuestos.

Usaremos a la función tsdiag( ) que proporciona gráficas para los residuales. En este caso lo estamos usando sobre el modelo que usa los datos originales.

plot(serie, xlim = c(1949,1962), col = "darkblue", lty = 1, lwd = 3,
ylab = "Pasajeros", xlab = "Tiempo en meses", legend = NULL,
main="Serie Original y Predicciones", ylim = c(100,700))
lines(ARMA11.pred$pred, col = "darkred", lty = 1, lwd = 2)
lines(ARMA11.pred$pred+1.96*ARMA11.pred$se, col = "darkred",
lty = 3, lwd = 2)
lines(ARMA11.pred$pred-1.96*ARMA11.pred$se, col = "darkred",
lty = 3, lwd = 2)
segments(1961-1/12, serie[length(serie)], 1961, ARMA11.pred$pred[1],
col = "darkred", lty = 1, lwd = 2)

tsdiag(ARMA11forma2)

Observamos que las autocorrelaciones de los residuos no se salen de las bandas de confianza lo cual sugiere un ruido blanco. Además todos los p-values del estadístico de Ljung Box están por arriba de un nivel de significancia de 0.05 hasta el lag 10 al menos, así que no rechazo la hipótesis nula de que las autocorrelaciones conjuntamente son cero, así que esto confirma que el residuo es de manera aproximada un ruido blanco. Finalmente se presentan los residuales y su histograma:

residuales<-ARMA11forma2$residuals
plot(residuales, xlab = "Tiempo", ylab = "Residuos", col = "darkred", lwd = 3)

hist (residuales, main = "Histograma residuales", prob=T, ylab = "Densidad",
breaks = 20, col="lightskyblue", border = "deepskyblue4", las = 1,
xlab = "Residuos")
grid(10, 10, lwd = 2)
curve(dnorm(x, mean = 0, sd = sd(residuales)), add = T, col = "darkblue",
lwd = 2)

# Una prueba formal de la normalidad es la de Lilliefors(Kolmogorov Sminrov):

lillie.test(residuales)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuales
## D = 0.084713, p-value = 0.01323

El p-value de la prueba de Liliefors es 0.01323 así que para un nivel de 0.05 rechazo la hipótesis de que los residuos se distribuyen normalmente. Realizaremos una prueba formal para probar si la varianza es constante:

ArchTest(residuales, lags=12)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  residuales
## Chi-squared = 20.167, df = 12, p-value = 0.06399

Para un nivel alfa de 0.05 no rechazo la hipótesis nula (usando s=12) y podemos concluir que que la variazna es constante.