Ajuste por regresión de una serie de tiempo con errores ARMA

Edwin Fernando Barrientos Arroyave
02/09/2021

Serie trimestral de los ingresos por acomodación, en hoteles licenciados que cuentan con 15 o más habitaciones, en Victoria (Australia), desde el trimestre 1 de 1998 al trimestre 3 de 2012 (en miles de dólares australianos)

plot of chunk unnamed-chunk-2

Serie trimestral de los ingresos por acomodación, en hoteles licenciados que cuentan con 15 o más habitaciones, en Victoria (Australia), desde el trimestre 1 de 1998 al trimestre 3 de 2012 (en miles de dólares australianos)

Modelo 1

Modelo aditivo de tendencia cuadrática y estacionalidad representada mediante funciones trigonométricas para las ondas sinusoidales armónicas en frecuencias angulares \( \pi \) y \( \frac{\pi}{2} \)

Ecuación

  • \( Y_{t}= \beta_0+\beta_{1}t+\beta_{2}t^2 + \sum_{j=1}^{2}[\alpha_{j}\sin(2 \pi F_{j}t) + \gamma_{j}\cos(2 \pi F_{j}t)] + E_{t}, \) con \( \{E_t\}_{t \in Z^{+} } un R.B \sim N(0,\sigma^{2}) \)

Análisis de residuales, test ACF y PACF

Residuales, ACF y PACF

plot of chunk unnamed-chunk-3

Ciclos y varianza constante

  • Existen ciclos en la serie de tiempo de los residuos, esto implica que \( \rho(k) =corr(E_t,E_{t+k}) \neq 0 \), para al menos un \( k \neq 0 \). También podría decirse que \( \rho(1) =corr(E_t,E_{t+1}) > 0 \) dando evidencia en contra del supuesto de independencia.
  • Por otro lado no hay evidencia en contra el supuesto de varianza costante, pues se observa un patrón aproximadamente costante a lo largo del tiempo \( Var[E_t]=\sigma^2_{t} \approx \sigma^2 \).
  • No hay evidencia en contra del supuesto de \( E[E_t]=0 \)

Análisis de residuales, test ACF y PACF

Residuales, ACF y PACF

plot of chunk unnamed-chunk-4

Test ACF y PACF

  • Para el test ACF con K = 1,2,..20 se va a probar respectivamente la siguiente hipótesis: \( H_0 :\rho(k) =corr(E_t,E_{t+k}) = 0 \) vs \( H_1 :\rho(k) \neq 0 \) y se tendrá en cuenta el estadístico \( \hat{\rho}(k) \stackrel{aprox}{\sim} N(0,\frac{1}{51}) \) y criterio de rechazo \( |\hat{\rho}(k)| \geq 2/\sqrt(51) \). Se detecta que \( \rho(k) \neq 0 \) para k=1,2,3 y 13

  • Por otro lado el test PACF con K = 1,2,..20 se va a probar respectivamente la siguiente hipótesis: \( H_0 :\phi_{kk} =corr(E_t,E_{t+k}) = 0 \) vs \( H_1 :\phi_{kk} \neq 0 \) y se tendrá en cuenta el estadístico \( \hat{\phi}_{kk} \stackrel{aprox}{\sim} N(0,\frac{1}{51}) \) y criterio de rechazo \( |\hat{\phi}_{kk}| \geq 2/\sqrt(51) \). Se detecta que \( \phi_{kk} \neq 0 \) para k=4

¿Et es estacionario en covarianza?

ACF muestral

plot of chunk unnamed-chunk-5

Estacionariedad en covarianza

No se cumple el supuesto de ruido blanco, sin embargo, para \( E_t \) no hay evidencia en contra de que su media es constante en cero y se puede asumir con varianza aprox. constante. Como el proceso observado es ergódigo según ACF, es decir \( \lim_{k \to \infty}\rho(k)=0 \) rápidamente, entonces se pude decir que el proceso estocástico para \( E_t \) del modelo global es estacionario en covarianza.

Test Durbin Watson

Postulación del Modelo

  • \[ Y_{t}= \beta_0+\beta_{1}t+\beta_{2}t^2 + \sum_{j=1}^{2}[\alpha_{j}\sin(2 \pi F_{j}t) + \gamma_{j}\cos(2 \pi F_{j}t)] + E_{t}, con \{E_t\}_{t \in Z^{+} } un R.B \sim N(0,\sigma^{2}) \]
  • donde \( E_{t}=\phi_1 E_{t-1} + a_t \), \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \) y \( |\phi_1|<1 \)
lag rho estimado Estadístico D-W VP rho>0 VP rho<0
1 0.6803441 0.6196832 0 1

Conclusión

entonces \( \rho(1) =corr(E_t,E_{t+1}) \) y se prueba que \( H_0 :\rho(1)=0 \) vs \( H_1 :\rho(1) > 0 \). Se elige contrastar la hipótesis nula contra la hipótesis \( H_1 :\rho(1) > 0 \), ya que el estadístico de la prueba \( d1 = 0.6196832 < 2 \), lo que sólo es posible si \( \rho(1) > 0 \). El valor P correspondiente, \( P(DW1 < 0.6803441) \approx 0.00 \) implica rechazo de \( H_0 \) y por tanto evidencia a favor de \( H1 \), es decir, que se detecta que en los residuales del modelo 1 hay autocorrelación positiva de orden 1, lo cual invalida supuesto de ruido blanco.

Identificación de modelos ARMA para los Et del modelo 1

ACF y PACF

plot of chunk unnamed-chunk-6

Con base en el análisis de la ACF y PACF muestrales presentados, se propone un AR(4), pues la ACF presenta un patron tipo cola y la PACF un patrón que puede ser tipo corte, tomando como último rechazo válido en k=4.

Modelo ARMA identificado

\( E_t=\phi_{1}E_{t-1}+\phi_{2}E_{t-2}+\phi_{3}E_{t-3}+\phi_{4}E_{t-4}+a_t \) con \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \)

Identificación de modelos ARMA para los Et del modelo 1

Función de R: armasubset 

Reordering variables and trying again:

plot of chunk unnamed-chunk-7

Con la función armasubsets() de la librería TSA se identifica en el primer renglón superior un AR(4)(solo con los parámetros \( \phi_1 \) y \( \phi_4 \))

Modelo ARMA identificado

\( E_t=\phi_{1}E_{t-1}+\phi_{4}E_{t-4}+a_t \) con \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \)

Identificación de modelos ARMA para los Et del modelo 1

Función de R: auto.arima()

alt text

Con la función auto.arima() de la librería forecast aplicada sobre el vector de valores de los residuales (sin convertirlos en objeto ts()) y usando criterio AIC, se identifica un ARMA(2,2)

Modelo ARMA identificado

\( E_t=\phi_{1}E_{t-1}+\phi_{2}E_{t-2}+a_t+ \theta_{1}a_{t-1}+\theta_{2}a_{t-2} \) con \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \)

Modelos de regresión tendencia, estacionalidad y errores estructurales ARMA

Modelo 2: Cuadrático – estacional con errores ARMA(2,2)

\( Y_{t}= \beta_0+\beta_{1}t+\beta_{2}t^2 + \sum_{j=1}^{2}[\alpha_{j}\sin(2 \pi F_{j}t) + \gamma_{j}\cos(2 \pi F_{j}t)] + E_{t}, \)

con \( E_t=\phi_{1}E_{t-1}+\phi_{2}E_{t-2}+\phi_{3}E_{t-3}+\phi_{4}E_{t-4}+a_t \), \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \)

Modelo 3: Cuadrático – estacional con errores siguiendo modelo ARMA(1,0)(1,0)[4]

\( Y_{t}= \beta_0+\beta_{1}t+\beta_{2}t^2 + \sum_{j=1}^{2}[\alpha_{j}\sin(2 \pi F_{j}t) + \gamma_{j}\cos(2 \pi F_{j}t)] + E_{t}, \)

con \( E_t=\Phi_{1}E_{t-4}+\phi_{1}E_{t-1}+\phi_{1}\Phi_{1}E_{t-5}+a_t \), \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \)

Modelo 4: Cuadrático – estacional con errores AR(4), sólo con los coeficientes autorregresivos \( \phi_1 \) y \( \phi_4 \)

\( Y_{t}= \beta_0+\beta_{1}t+\beta_{2}t^2 + \sum_{j=1}^{2}[\alpha_{j}\sin(2 \pi F_{j}t) + \gamma_{j}\cos(2 \pi F_{j}t)] + E_{t}, \)

con \( E_t=\phi_{1}E_{t-1}+\phi_{4}E_{t-4}+a_t \), \( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \)

Parámetros estimados en modelos con erroes ARMA

Modelo 1

Estimación s.e Valor To P(t>|To|)
beta_o 79730.59841 2963.22628 26.90669 0.00000
beta_1 1750.84660 262.89587 6.65985 0.00000
beta_2 5.73786 4.90169 1.17059 0.24852
sen1 6256.17291 1327.49717 4.71276 0.00003
sen2 8019.12167 1355.73986 5.91494 0.00000
sen3 -3463.68370 948.28236 -3.65259 0.00073

Modelo 3

Estimación s.e Valor To P(t>|To|)
ar1 0.69580 0.09974 6.97592 0.00000
sar1 -0.10456 0.14850 -0.70407 0.48518
beta_o 79166.76531 5045.77892 15.68970 0.00000
beta_1 1750.10736 446.20718 3.92219 0.00031
beta_2 6.44905 8.28632 0.77828 0.44067
sen1 6419.64088 685.60117 9.36352 0.00000
sen2 8058.18107 687.54295 11.72026 0.00000
sen3 -3458.65243 349.91002 -9.88441 0.00000

Modelo 2

Estimación s.e Valor To P(t>|To|)
ar1 1.50448 0.20467 7.35072 0.00000
ar2 -0.74089 0.16768 -4.41838 0.00007
ma1 -0.97269 0.19495 -4.98933 0.00001
ma2 0.55320 0.17165 3.22284 0.00249
beta_o 79824.57609 4418.75212 18.06496 0.00000
beta_1 1661.82930 397.55099 4.18017 0.00015
beta_2 8.47511 7.48751 1.13190 0.26425
sen1 6435.72464 592.88441 10.85494 0.00000
sen2 8082.03076 602.19419 13.42097 0.00000
sen3 -3420.39479 465.64257 -7.34554 0.00000

Modelo 4

Estimación s.e Valor To P(t>|To|)
ar1 0.75788 0.10420 7.27335 0.00000
ar4 -0.19842 0.10406 -1.90683 0.06324
beta_o 79725.99675 4420.92801 18.03377 0.00000
beta_1 1676.44220 394.95828 4.24461 0.00011
beta_2 8.11559 7.38923 1.09830 0.27819
sen1 6450.16315 629.52412 10.24609 0.00000
sen2 8053.33038 633.17495 12.71897 0.00000
sen3 -3444.06770 323.66520 -10.64083 0.00000

Gráficos de residuales(Modelos: 2,3 y 4)

Residuos vs tiempo

plot of chunk unnamed-chunk-9

En todas las gráficas de residuales de ajuste de los modelos 2 a 4 correspondiente a \( \hat{a}_t \), no se encuentra evidencia en contra \( E[a_t]=0 \) y \( Var[a_t]=\sigma^2_{a} \), \( \forall_t \).

Residuos vs valores ajustados

plot of chunk unnamed-chunk-10

Gráficos ACF y PACF(Modelos: 2,3 y 4)

ACF´s de los residuos de ajuste

plot of chunk unnamed-chunk-11 En los test ACF, donde se define \( \rho(k)=corr(a_t,a_{t+k}) \) con \( k=1,2,..,20 \) no se halla evidencia en contra del supuesto R.B.

PACF´s de los residuos de ajuste

plot of chunk unnamed-chunk-12 En los test PACF, donde se define \( \phi_{kk}=corr(a_t,a_{t+k}|a_{t+1},...,a_{t+k-1}) \) con \( k=1,2,..,20 \) no se halla evidencia en contra del supuesto R.B.

Gráficos de normalidad(Modelos: 2,3 y 4)

Gráficos de probabilidad normal

plot of chunk unnamed-chunk-13

Test a probar

\( H_0:a_t \sim N(0,\sigma^2_a) \) vs \( H_1:a_t \nsim N(0,\sigma^2_a) \)

Estadístico W Valor P
modelo2 0.9697316 0.2155589
modelo3 0.9707561 0.2376867
modelo4 0.9640358 0.1239596

En los modelos 2,3 y 4 se observan fuestes desviaciones a lo largo de toda la nubes de puntos.

Comparación de los ajustes

Comparación de gráficos de ajuste

plot of chunk unnamed-chunk-14

Comparación de medidas de ajuste

Modelo1 Modelo2 Modelo3 Modelo4
AIC 51032490 26783180 28763331 26999194
BIC 64054615 39117245 38944315 36555750

Gráficamente, se observa una mejora en los ajustes al incluir las estructuras ARMA en los errores estructurales \( E_t \), pero es difícil distinguir por gráficos en cuál de los modelos 2 a 4 es mejor el ajuste. Con los criterios de información AIC y BIC se encuentra que por ambos criterios el mejor ajuste es con el modelo 4

Pronósticos ex-post

Comparación precisión pronósticos puntuales y por I.P

RMSE MAE MAPE Amp Media I.P Cobertura
mod1 9226.321 8416.517 4.093124 32149.49 1
mod2 4689.718 4078.036 1.977711 24576.38 1
mod3 6911.309 6167.917 2.992077 25282.03 1
mod4 5004.117 4261.396 2.056165 25897.02 1

Comparación precisión pronósticos puntuales y por I.P

Considerando que en los modelos 2,3 y 4 no se imcumple el supuesto de (\( \{a_t\}_{t\in z^+} \) un R.B \( \sim N(0,\sigma^2_{a}) \) ).Los I.P de los modelos logran una cobertura de \( 100\% \) de los valores reales. EL modelo 2 tiene un I.P mas estrecho por otro lado, usando las medidas de MAE, MAPE y RMSE de pronóstico se halla que el error promedio en pronóstico es un poco menor en el modelo 2, seguido por el modelo 4.

Pronósticos ex-post

Comparación gráfica de pronósticos

plot of chunk unnamed-chunk-15

Terminar

  • En terminos generales los modelos están subestimando la mayoría de tiempo los pronósticos reales, a excepción de los modelos 2 y 4 que en algunos periodos sobre estiman dichos pronósticos. También se observa en la gráfica que el modelo 2 es el más cercano a la realidad.

  • El modelo que mejor pronóstica es el modelo 2 ya que:

    • Tiene I.P mas estrecho que el resto de modelos
    • Tiene un error promedio en pronóstico menor que el resto

Modelo a recomendar

Teniendo como primer criterio de selección los resultados en el análisis de residuales, en segundo lugar los resultados de pronósticos y por último los de ajuste, se concluye que el mejor modelo por el momento, es el modelo 2.