Análisis de una serie temporal. Tráfico comercial en el aeropuerto de Valencia.

En este caso práctico analizamos la serie temporal correspondiente al tráfico comercial en el aeropuerto de Valencia. Cargamos los dos paquetes de R que vamos a utilizar en el análisis.

library(readr)
library(forecast)

Ejemplo: Tráfico comercial, aeropuerto de Valencia. Fuente: Ministerio de transportes, movilidad y agenda urbana (https://www.fomento.gob.es/BE/?nivel=2&orden=03000000).

Dat_Pasajeros <- read.csv(file="Pasajeros.csv",header=TRUE,sep=";")
attach(Dat_Pasajeros)
T <- length(Total)
Pasajeros_ts <- ts(Total[121:T],start=c(2010,1),end=c(2019,12),frequency=12)
# Creamos una serie temporal con las observaciones a partir de enero de 2010. 
# Como se trata de datos mensuales, definimos frequency=12

1 Descripción gráfica de la serie temporal

Es habitual comenzar el análisis de una serie con la representación gráfica de los valores observados de la variable de interés en función del tiempo:

plot(Pasajeros_ts,ylab="Num Pasajeros")

A partir del gráfico temporal podemos apreciar una cierta evolución en el largo plazo (tendencia): durante los primeros cuatro años se observa una tendencia decreciente. A partir de 2014, la serie temporal toma valores cada vez mayores, es decir, la serie presenta una tendencia creciente. Por otro lado, se observa un comportamiento cíclico que se repite año tras año (estacionalidad), con un número mayor de pasajeros durante los meses de verano. La longitud del ciclo estacional es \(c = 12\).

En este ejemplo, la estacionalidad de la serie se observa claramente en el gráfico temporal. No obstante, el diagrama de cajas por mes nos permite también valorar la presencia de estacionalidad.

Mes_ord <- factor(Mes[121:T],levels = c("Enero","Febrero","Marzo","Abril",
"Mayo","Junio","Julio","Agosto","Septiembre","Octubre","Noviembre","Diciembre"))
# Ordenamos los meses para que los represente por orden temporal

boxplot(Total[121:T] ~ Mes_ord, ylab="", xlab="")

2 Análisis de la serie mediante suavizado exponencial

Dadas las características de la serie temporal: tendencia y estacionalidad, el método adecuado para su análisis es el método de Holt-Winters.

Vamos a empezar analizando la serie con el método de Holt-Winters aditivo pues, a partir del gráfico temporal, podríamos asumir que el efecto de la estacionalidad es aditivo (no parece aumentar con el nivel). No obstante, analizaremos también la serie con el método de Holt-Winters multiplicativo y la serie transformada (utilizando la transformación logarítmica) con Holt-Winters aditivo.

Utilizamos los datos hasta diciembre de 2018 para el ajuste y reservamos las observaciones de 2019 para valorar la capacidad predictiva del método seleccionado. La predicción para el año 2019 la realizaremos utilizando el método que nos proporciones un mejor ajuste.

2.1 Holt-Winters aditivo

insample <- window(Pasajeros_ts,start=c(2010,1),end=c(2018,12))  
      # ajuste desde enero de 2010 hasta diciembre de 2018
outsample <- window(Pasajeros_ts,start=c(2019,1),end=c(2019,12)) 
      # reservamos 2019 para valorar la predicción

fitPasajeros <- HoltWinters(insample,seasonal="additive")
fitPasajeros$coefficients

##          a          b         s1         s2         s3         s4         s5 
##  664886.05    6144.59 -137983.00 -142780.65  -25444.56   38806.44   53157.14 
##         s6         s7         s8         s9        s10        s11        s12 
##   67880.47  164551.61  194081.05   90136.64   64831.64  -77773.59 -100537.05

fitPasajeros$alpha

##    alpha 
## 0.325429

fitPasajeros$beta

##       beta 
## 0.06326808

fitPasajeros$gamma

## gamma 
##     1

fitval <- fitted(fitPasajeros)  
# fitval contiene la serie de valores ajustados en la primera columna (fitval[,1] = xhat)

plot(fitPasajeros,ylab="Num Pasajeros",main="Ajuste HW aditivo")

# Valoramos la bondad del ajuste

insamplecut <- window(insample,start=c(2011,1),end=c(2018,12))
# El año 2010 se utiliza para calcular las condiciones iniciales. 
# El ajuste pues se obtiene a partir de enero de 2011.

rmse <- sqrt(mean((insamplecut-fitval[,1])^2))
mape <- 100*mean(abs(insamplecut-fitval[,1])/insamplecut)
rmse; mape

## [1] 18481.51

## [1] 3.480071

2.2 Holt-Winters multiplicativo

fitPasajeros_mult <- HoltWinters(insample,seasonal="multiplicative")
fitPasajeros_mult$coefficients

##            a            b           s1           s2           s3           s4 
## 7.252817e+05 3.680679e+03 7.355365e-01 7.549979e-01 9.829124e-01 1.126145e+00 
##           s5           s6           s7           s8           s9          s10 
## 1.108105e+00 1.105080e+00 1.252713e+00 1.319936e+00 1.118115e+00 1.044445e+00 
##          s11          s12 
## 8.020068e-01 7.781101e-01

fitPasajeros_mult$alpha

##     alpha 
## 0.8073643

fitPasajeros_mult$beta

##       beta 
## 0.01386687

fitPasajeros_mult$gamma

## gamma 
##     1

fitval_mult <- fitted(fitPasajeros_mult)  

plot(fitPasajeros_mult,ylab="Num Pasajeros",main="Ajuste HW multiplicativo")

# Valoramos la bondad del ajuste

rmse_mult <- sqrt(mean((insamplecut-fitval_mult[,1])^2))
mape_mult <- 100*mean(abs(insamplecut-fitval_mult[,1])/insamplecut)
rmse_mult; mape_mult

## [1] 23384

## [1] 3.852084

2.3 Holt-Winters aditivo aplicado a la serie transformada

loginsample <- log(insample) 

fitPasajeros_log <- HoltWinters(loginsample,seasonal="additive")
fitPasajeros_log$coefficients

##            a            b           s1           s2           s3           s4 
## 13.378855434  0.008592261 -0.193608764 -0.218809475  0.003296337  0.094101281 
##           s5           s6           s7           s8           s9          s10 
##  0.101634949  0.107006599  0.233074551  0.267322754  0.134865747  0.108210330 
##          s11          s12 
## -0.105720853 -0.135427300

fitPasajeros_log$alpha

##      alpha 
## 0.07440285

fitPasajeros_log$beta

##      beta 
## 0.2434083

fitPasajeros_log$gamma

## gamma 
##     1

fitval_log <- fitted(fitPasajeros_log)  

plot(fitPasajeros_log,ylab="Log(Num Pasajeros)",main="Ajuste HW aditivo a la serie de los logaritmos")

# Valoramos la bondad del ajuste. Para ello, volvemos previamente a la escala original

fitval_ori <- exp(fitval_log[,1])

rmse_log <- sqrt(mean((insamplecut-fitval_ori)^2))
mape_log <- 100*mean(abs(insamplecut-fitval_ori)/insamplecut)
rmse_log; mape_log

## [1] 25480.46

## [1] 3.991273

El método con menor error de ajuste (tanto RMSE como MAPE) es Holt-Winters con estacionalidad aditiva. Este será, por tanto, el método utilizado para calcular la predicción para el año 2019.

pred <- predict(fitPasajeros,12)
# pred contiene las predicciones puntuales para los 12 meses de 2019
ts.plot(insample,pred,lty=1:2)

# Valoramos la capacidad predictiva del método

rmse_pred <- sqrt(mean((outsample-pred)^2))
mape_pred <- 100*mean(abs(outsample-pred)/outsample)
rmse_pred;mape_pred

## [1] 22553.38

## [1] 2.832105

Podemos también representar gráficamente los valores reales de 2019 que habíamos reservado junto con la predicción puntual:

plot(pred, col="red",xaxt="n",xlab="Año 2019")
points(outsample,pch=19)

y calcular el intervalo de predicción al 95%:

pred <- predict(fitPasajeros,n.ahead=12,prediction.interval=TRUE,level=0.95) 
plot(fitPasajeros, pred)

3 Análisis de la serie mediante la metodología Box-Jenkins

Dadas las características de la serie temporal: tendencia y estacionalidad, el primer paso del análisis es determinar la transformación estacionaria de la serie.

Calculamos un diferencia estacional (\(D = 1\)):

d12insample <- diff(insample,12)
plot(d12insample)

Parece que hemos quitado la estacionalidad, pero todavía queda la tendencia. Calculamos pues una diferencia regular (\(d = 1\)):

dd12insample <- diff(d12insample)
plot(dd12insample)

Podemos asumir que la serie diferenciada con \(d = 1\) y \(D = 1\) ya es estacionaria. Pasamos a examinar el correlograma y el correlograma parcial:

acf(dd12insample,lag.max=50)

pacf(dd12insample,lag.max=50)

Si nos fijamos en los primeros retardos, podemos pensar:

• La función de autocorrelación tiene el primer coeficiente significativo, mientras que la función de autocorrelación parcial muestra decrecimiento: (p,d,q) = (0,1,1)

• La función de autocorrelación decrece y la función de autocorrelación parcial tiene el primer coeficiente significativo: (p,d,q) = (1,1,0)

• Las dos funciones muestran decrecimiento a partir del primer coeficiente: (p,d,q) = (1,1,1)

Si nos fijamos en los retardos estacionales (Lag = 1, 2, 3, 4 ciclos estacionales), podemos pensar:

• No hay ningún coeficiente significativo: (P,D,Q) = (0,1,0)

• La función de autocorrelación tiene el primer coeficiente significativo, mientras que la función de autocorrelación parcial muestra decrecimiento: (P,D,Q) = (0,1,1)

Veamos el ajuste proporcionado por los distintos modelos:

Pasajeros_model1 <- arima(insample, order=c(0,1,1), seasonal=list(order=c(0,1,0), period=12))
Pasajeros_model1

## 
## Call:
## arima(x = insample, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0), 
##     period = 12))
## 
## Coefficients:
##           ma1
##       -0.4055
## s.e.   0.0965
## 
## sigma^2 estimated as 269599887:  log likelihood = -1056.98,  aic = 2117.96

Pasajeros_model2 <- arima(insample, order=c(0,1,1), seasonal=list(order=c(0,1,1), period=12))
Pasajeros_model2

## 
## Call:
## arima(x = insample, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), 
##     period = 12))
## 
## Coefficients:
##           ma1     sma1
##       -0.3712  -0.2461
## s.e.   0.0980   0.1275
## 
## sigma^2 estimated as 258053809:  log likelihood = -1055.26,  aic = 2116.52

Pasajeros_model3 <- arima(insample, order=c(1,1,0), seasonal=list(order=c(0,1,0), period=12))
Pasajeros_model3

## 
## Call:
## arima(x = insample, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 0), 
##     period = 12))
## 
## Coefficients:
##           ar1
##       -0.3601
## s.e.   0.0955
## 
## sigma^2 estimated as 273628695:  log likelihood = -1057.66,  aic = 2119.33

Pasajeros_model4 <- arima(insample, order=c(1,1,0), seasonal=list(order=c(0,1,1), period=12))
Pasajeros_model4

## 
## Call:
## arima(x = insample, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 1), 
##     period = 12))
## 
## Coefficients:
##           ar1     sma1
##       -0.3353  -0.2627
## s.e.   0.0972   0.1303
## 
## sigma^2 estimated as 260716590:  log likelihood = -1055.79,  aic = 2117.58

Pasajeros_model5 <- arima(insample, order=c(1,1,1), seasonal=list(order=c(0,1,0), period=12))
Pasajeros_model5

## 
## Call:
## arima(x = insample, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), 
##     period = 12))
## 
## Coefficients:
##           ar1      ma1
##       -0.0166  -0.3908
## s.e.   0.3056   0.2904
## 
## sigma^2 estimated as 269590385:  log likelihood = -1056.98,  aic = 2119.96

Pasajeros_model6 <- arima(insample, order=c(1,1,1), seasonal=list(order=c(0,1,1), period=12))
Pasajeros_model6

## 
## Call:
## arima(x = insample, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), 
##     period = 12))
## 
## Coefficients:
##          ar1      ma1     sma1
##       -0.038  -0.3369  -0.2474
## s.e.   0.313   0.3020   0.1281
## 
## sigma^2 estimated as 2.58e+08:  log likelihood = -1055.25,  aic = 2118.51

Pasajeros_model <- auto.arima(insample)
Pasajeros_model

## Series: insample 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.3712  -0.2461
## s.e.   0.0980   0.1275
## 
## sigma^2 estimated as 263622661:  log likelihood=-1055.26
## AIC=2116.52   AICc=2116.78   BIC=2124.18

accuracy(Pasajeros_model)

##                    ME    RMSE      MAE          MPE     MAPE      MASE
## Training set 815.3097 15066.8 10530.45 -0.003326076 2.459859 0.2577526
##                      ACF1
## Training set -0.005042876

El modelo de menor AIC es (p,d,q)(P,D,Q) = (0,1,1)(0,1,1), que coincide con el modelo proporcionado por la función auto.arima. El MAPE asociado a este modelo es 2.4599. La ecuación del modelo es: \[ \bigtriangledown \bigtriangledown_{12} x_t = (1 - 0.3712 B) (1 - 0.2461 B^{12})\epsilon_t \] \[ (1 - B) (1 - B^{12}) x_t = (1 - 0.3712 B) (1 - 0.2461 B^{12})\epsilon_t \]

Veamos a continuación la representación gráfica del ajuste obtenido. Línea negra: valores reales, línea roja: valores ajustados.

fitval <- Pasajeros_model$fitted # Valores ajustados

plot(insample,ylab="Num Pasajeros")
lines(fitval,col="red")

Antes de pasar a la predicción, comprobamos que el modelo es válido. Como muestran las siguientes salidas, los residuos del modelo pueden considerarse ruido blanco.

checkresiduals(Pasajeros_model,plot=TRUE)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 10.68, df = 20, p-value = 0.9541
## 
## Model df: 2.   Total lags used: 22

La predicción obtenida para los 12 meses de 2019 junto con el error de predicción vienen dados por:

pred <- forecast(Pasajeros_model,h=12)$mean
pred # Predicción puntual

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2019 539654.2 539193.2 660745.2 733403.9 752033.1 768074.1 865913.7 900445.2
##           Sep      Oct      Nov      Dec
## 2019 802007.1 785309.0 653201.9 639994.1

plot(forecast(Pasajeros_model,h=12))

rmse_pred <- sqrt(mean((outsample-pred)^2))
mape_pred <- 100*mean(abs(outsample-pred)/outsample)
rmse_pred;mape_pred

## [1] 20479.29

## [1] 2.486239

Finalmente, representamos gráficamente los valores reales de 2019 que habíamos reservado junto con la predicción puntual:

plot(pred, col="red",xaxt="n",xlab="Año 2019")
points(outsample,pch=19)

Si comparamos ambas metodologías, vemos que el error de ajuste correspondiente al modelo sARIMA es menor que el obtenido con Holt-Winters aditivo y, por tanto, como predicción para el año 2019 deberíamos haber tomado las obtenidas con la metodología Box-Jenkins. Además, como habíamos reservado las observaciones de 2019 para valorar la capacidad predictiva del modelo, comprobamos que las predicciones obtenidas con el modelo sARIMA son, también, más precisas (menor error de predicción).