Taller SARIMA

Author

Daniela Alcázar

Librerías usadas

library(fpp2) 
library(forecast)
library(ggplot2)
library(tseries)
library(Metrics)
library(car)
library(nortest)
library(lmtest)

Exploración de la serie

Vamos a considerar la serie a10, la cual contiene los registros mensuales del gasto gubernamental (en millones de dólares) en subsidios para medicamentos antidiabéticos en Australia, a partir de julio de 1991 hasta junio de 2008.

Para empezar vamos a leer y a renombrar la serie. Esta serie hace parte de la librería fpp2, la cuál se encuentra disponible en la red CRAN.

if (!requireNamespace("fpp2", quietly = TRUE)) {
    install.packages("fpp2")}
library(fpp2)

data <- a10

Ahora identificaremos las características de la serie.

typeof(data)

[1] "double"

Los registros de la serie son de tipo decimal.

length(data)

[1] 204

La serie contiene 204 registros.

frequency(data)

[1] 12

La frecuencia de la serie es mensual.

start(data)

[1] 1991    7

El primer registro que se tiene en la serie es de julio del 1991.

end(data)

[1] 2008    6

El último registro de la serie es de junio del 2008.

Identificación de componentes

autoplot(data)

A partir de la gráfica de la serie podemos llegar a las siguientes conjeturas:

La serie al inicio de los registros toma valores cercanos al 5, con el tiempo se va evidenciando un paulatino crecimiento conforme pasan los años, para los últimos años la serie toma valores alrededor del 20. La media de la serie no es estacionaria, es decir, esta no permanece constante y parece en cambio depender del tiempo.
La serie parece tener un componente estacional, ya que año a año se repite un comportamiento característico en el cual los valores crecen hasta alcanzar un valor pico a final de cada año que luego baja dramáticamente para volver a empezar su patrón de comportamiento. Como la longitud de este comportamiento es fijo podemos pensar que es un comportamiento estacional y no cíclico.
Con el pasar del tiempo los valores parecen cambiar con más fuerza entre un periodo y otro, la variación entre datos de la serie es menor al inicio de los registros y mayor al final. Esto nos dice que la varianza de la serie no es estacionaria, es decir, que parece depender del tiempo.

Recurrimos al análisis de la gráfica de la función de autocorrelación ACF y de la función de autocorrelaciones parciales PACF para corroborar nuestras hipótesis.

ggAcf(data, lag.max = 50)
ggPacf(data, lag.max = 50)

La ACF tiene valores altos en todos los rezagos, los rezagos tienen forma de U con los valores más altos en los periodos estacionales (rezagos múltiplos de 12). La PACF tiene un comportamiento más irregular, con un valor fuera de las bandas en el primer rezago estacional (rezago 12)

Para concluir, la serie no es estacionaria. Depende del tiempo, pues, además del componente estacional, su media y su varianza no son estacionarias, ya que varían en función del tiempo.

Otras herramientas gráficas

Existen funciones como decompose y ggseasonalplot que pueden ser usadas para identificar componentes presentes en una serie de tiempo.

A continuación se presentan algunas formas de usarlas

plot(decompose(data, type = 'additive')) # Series con varianza estacionaria
plot(decompose(data, type = 'multiplicative')) # Series con varianza no estacionaria
ggseasonplot(data, polar=TRUE) +
    ylab("Millones $")
ggseasonplot(data, polar=FALSE) +
    ylab("Millones $")

Transformación y diferenciación

Para poder pensar en un modelo SARIMA debemos garantizar que la serie es estacionaria.

Primero se estabiliza la varianza de la serie por medio de una transformación, en este caso vamos a aplicar la transformación Box-Cox.

lambda <- BoxCox.lambda(data)
lambda

[1] 0.1313326

En este caso identificamos un valor de 0.1313326 para lambda, el parámetro de la transformación Box-Cox.

data_T <- BoxCox(data, lambda)
plot(data_T)

Es evidente que la transformación ayuda a alcanzar la estacionariedad en varianza, por esta razón decidimos trabajar en adelante con la serie transformada.

Continuamos con la diferenciación de la serie. Para series de tiempo que necesiten de diferenciación estacional y regular (series con tendencia y estacionalidad), se pueden aplicar en el orden que el analista lo decida. En este caso se decide aplicar una diferenciación estacional primero.

Para empezar, identificamos el número de diferenciaciones estacionales que debemos aplicar.

nsdiffs(data_T)

[1] 1

Ahora diferenciamos la serie con respecto a la longitud del rezago estacional.

data_T_DS <- diff(data_T, lag = frequency(data))
plot(data_T_DS)

Con la varianza estabilizada y sin el componente estacional procedemos a estabilizar la media. Antes, para confirmar que esta no es estacionaria aplicamos la prueba de la raíz unitaria Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

Recordemos que la prueba KPSS presenta las siguientes hipótesis:

\[ H_0 = La\;serie\;es\;estacionaria \] \[ H_1 = La\;serie\;no\;es\;estacionaria \]

kpss.test(data_T_DS)

Warning in kpss.test(data_T_DS): p-value greater than printed p-value


    KPSS Test for Level Stationarity

data:  data_T_DS
KPSS Level = 0.14853, Truncation lag parameter = 4, p-value = 0.1

El valor p de la prueba KPSS tiene un valor mayor que 0.1, por lo tanto, aceptamos la hipótesis nula que dice que la serie es estacionaria. Con la diferenciación estacional también se eliminó la tendencia, como la media ya es estacionaria no hay necesidad de aplicar una diferenciación regular. Por medio de la diferenciación estacional se observa un proceso estacionario a partir del cuál podemos empezar a concluir sobre los modelos SARIMA.

Identificación de parámetros

Como estamos trabajando con una serie con componente estacional, debemos pensar en los parámetros p, q, P y Q. Para esto gracficamos la ACF y la PACF de la serie estacionaria.

ggAcf(data_T_DS, lag.max = 50)
ggPacf(data_T_DS, lag.max = 50)

Para determinar los parámetros P y Q analizamos los rezagos estacionales, estos son los múltiplos de la frecuencia de la serie, es decir 12, 24, 36 y 48.
Para ambos casos, el ACF y el PACF, solo los dos primeros rezagos estacionales (12 y 24) se salen de las bandas. El ACF es el que más parece tener un comportamiento especial, pues se podría decir que hay un posible comportamiento sinusoidal. Esto nos lleva a considerar un posible ARMA(2, 2) estacional O AR(2) estacional.

Continuando con los patrones regulares, estudiamos los rezagos que hacen parte de un periodo estacional. De estos se puede notar que en la ACF y en la PACF son los tres primeros rezagos regulares los que están por fuera de las bandas. También se podría pensar que el ACF decae lentamente. En consecuencia, podemos pensar en los modelos ARMA(3, 3) y AR(3).

Identificación modelos

De acuerdo al procedimiento realizado y con base en el análisis previamente hecho, se identifican los siguientes modelos:

Modelo 1: SARIMA(3,0,3)(2,1,2)
Modelo 2: SARIMA(3,0,3)(2,1,0)

Modelo 3: SARIMA(3,0,0)(2,1,2)
Modelo 4: SARIMA(3,0,0)(2,1,0)

Ahora, por último, usaremos la función auto.arima para identificar algún modelo que no hayamos podido identificar.

auto.arima(data_T, approximation = FALSE)

Series: data_T 
ARIMA(3,0,0)(2,1,1)[12] with drift 

Coefficients:
         ar1     ar2     ar3    sar1     sar2     sma1   drift
      0.0008  0.2681  0.3164  0.0645  -0.2576  -0.7389  0.0126
s.e.  0.0720  0.0678  0.0754  0.1187   0.1019   0.1041  0.0003

sigma^2 = 0.006082:  log likelihood = 214.12
AIC=-412.23   AICc=-411.44   BIC=-386.17

De esto identificamos nuestro último modelo, el modelo 5:

Modelo 5: SARIMA(3,0,0)(2,1,1)

Ajuste de los modelos

Vamos a estimar los coeficientes para cada modelo y a probar su significancia usando la prueba de la normal.

Recordemos que esta prueba presenta las siguientes hipótesis:

Para el caso de los coeficientes $\phi_i$

\[ H_0: \phi_i=0 \]

\[ H_1: \phi_i\neq 0 \]

Para el caso de los coeficientes $\theta_j$

\[ H_0: \theta_j=0 \]

\[ H_1: \theta_j\neq 0 \]

Para los coeficientes estacionales es igual, la diferencia es que en el subíndice del órden también se agrega una “s”

Modelo 1: SARIMA(3,0,3)(2,1,2)

modelo1 <- Arima(data_T_DS, order = c(3,0,3), seasonal = c(2,0,2))
summary(modelo1)

Series: data_T_DS 
ARIMA(3,0,3)(2,0,2)[12] with non-zero mean 

Coefficients:
          ar1     ar2     ar3    ma1     ma2      ma3     sar1     sar2
      -0.0964  0.0091  0.7184  0.117  0.2825  -0.4280  -0.3193  -0.1826
s.e.   0.1291  0.1204  0.1199  0.155  0.1303   0.1661   0.3858   0.1467
         sma1     sma2    mean
      -0.3863  -0.3268  0.1511
s.e.   0.3901   0.3613  0.0037

sigma^2 = 0.006031:  log likelihood = 216.85
AIC=-409.71   AICc=-407.96   BIC=-370.62

Training set error measures:
                      ME       RMSE        MAE       MPE     MAPE      MASE
Training set 0.001646662 0.07540344 0.05725127 -25.13397 74.37486 0.5004162
                     ACF1
Training set -0.004962855

matriz <- cbind(modelo1$coef, sqrt(diag(modelo1$var.coef)))
Z = matriz[,1]/matriz[,2]

valor_p <- 2*pnorm(-abs(Z))
valor_p

         ar1          ar2          ar3          ma1          ma2          ma3 
4.554584e-01 9.394704e-01 2.079432e-09 4.502219e-01 3.011454e-02 9.985373e-03 
        sar1         sar2         sma1         sma2    intercept 
4.079341e-01 2.131559e-01 3.219872e-01 3.657055e-01 0.000000e+00

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\theta_3$ (ma3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

$\phi_{2s}$ (sar2) ⟶ Como el valor p es mayor que la significancia $\alpha = 0.05$, se concluye que no es significativo.
$\theta_{2s}$ (sma2) ⟶ Como el valor p es mayor que la significancia $\alpha = 0.05$, se concluye que no es significativo.

De este modelo se concluye que $\phi_{2s}$ y $\theta_{2s}$ no son significativos. Se propone el modelo 1.1: SARIMA(3,0,3)(1,1,1)

Modelo 1.1: SARIMA(3,0,3)(1,1,1)

modelo1_1 <- Arima(data_T_DS, order = c(3,0,3), seasonal = c(1,0,1))
summary(modelo1_1)

Series: data_T_DS 
ARIMA(3,0,3)(1,0,1)[12] with non-zero mean 

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3    sar1     sma1    mean
      -0.0906  0.0004  0.7332  0.1148  0.2874  -0.4543  0.1425  -0.9039  0.1516
s.e.   0.1237  0.1133  0.1054  0.1487  0.1251   0.1416  0.1127   0.1070  0.0034

sigma^2 = 0.006042:  log likelihood = 214.35
AIC=-408.7   AICc=-407.48   BIC=-376.12

Training set error measures:
                      ME       RMSE        MAE       MPE     MAPE      MASE
Training set 0.001254799 0.07588509 0.05774877 -26.51159 75.61896 0.5047647
                     ACF1
Training set -0.008222121

matriz <- cbind(modelo1_1$coef, sqrt(diag(modelo1_1$var.coef)))
Z = matriz[,1]/matriz[,2]
valor_p <- 2*pnorm(-abs(Z))
valor_p

         ar1          ar2          ar3          ma1          ma2          ma3 
4.640190e-01 9.972093e-01 3.445067e-12 4.400878e-01 2.157917e-02 1.330904e-03 
        sar1         sma1    intercept 
2.060934e-01 3.037231e-17 0.000000e+00

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\theta_3$ (ma3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

$\phi_{1s}$ (sar1) ⟶ Como el valor p es mayor que la significancia $\alpha = 0.05$, se concluye que no es significativo.
$\theta_{1s}$ (sma1) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

De este modelo se concluye que $\theta_{1s}$ no es significativo. Se propone el modelo 1.2: SARIMA(3,0,3)(1,1,0)

Modelo 1.2: SARIMA(3,0,3)(1,1,0)

modelo1_2 <- Arima(data_T_DS, order = c(3,0,3), seasonal = c(1,0,0))
summary(modelo1_2)

Series: data_T_DS 
ARIMA(3,0,3)(1,0,0)[12] with non-zero mean 

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3     sar1    mean
      -0.1431  0.0808  0.7031  0.1586  0.1892  -0.4372  -0.3831  0.1574
s.e.   0.1927  0.1057  0.2044  0.2474  0.1147   0.2074   0.0806  0.0116

sigma^2 = 0.008135:  log likelihood = 192.38
AIC=-366.77   AICc=-365.78   BIC=-337.45

Training set error measures:
                      ME       RMSE        MAE       MPE     MAPE      MASE
Training set 0.000349926 0.08829696 0.06770023 -34.23164 90.01941 0.5917474
                    ACF1
Training set -0.01249447

matriz <- cbind(modelo1_2$coef, sqrt(diag(modelo1_2$var.coef)))
Z = matriz[,1]/matriz[,2]  
valor_p <- 2*pnorm(-abs(Z)) 
valor_p

         ar1          ar2          ar3          ma1          ma2          ma3 
4.577163e-01 4.447907e-01 5.808117e-04 5.215047e-01 9.913401e-02 3.499163e-02 
        sar1    intercept 
1.986007e-06 7.607330e-42

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\theta_3$ (ma3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

$\phi_{1s}$ (sar1) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

De este modelo se concluye que todos los parámetros son significativos.
Modelo 2: SARIMA(3,0,3)(2,1,0)

modelo2 <- Arima(data_T_DS, order = c(3,0,3), seasonal = c(2,0,0))
summary(modelo2)

Series: data_T_DS 
ARIMA(3,0,3)(2,0,0)[12] with non-zero mean 

Coefficients:
          ar1     ar2     ar3     ma1     ma2      ma3     sar1     sar2
      -0.0839  0.0846  0.7053  0.0476  0.2035  -0.4009  -0.5439  -0.4771
s.e.   0.1343  0.1044  0.1213  0.1664  0.1223   0.1452   0.0804   0.0768
        mean
      0.1546
s.e.  0.0085

sigma^2 = 0.006732:  log likelihood = 208.16
AIC=-396.33   AICc=-395.11   BIC=-363.75

Training set error measures:
                       ME       RMSE        MAE       MPE     MAPE      MASE
Training set 0.0009854997 0.08010045 0.06216446 -31.16602 82.10343 0.5433609
                     ACF1
Training set 0.0007931419

matriz <- cbind(modelo2$coef, sqrt(diag(modelo2$var.coef)))
Z = matriz[,1]/matriz[,2]
valor_p <- 2*pnorm(-abs(Z))
valor_p

         ar1          ar2          ar3          ma1          ma2          ma3 
5.322621e-01 4.180915e-01 6.135138e-09 7.746339e-01 9.605035e-02 5.753360e-03 
        sar1         sar2    intercept 
1.330172e-11 5.243850e-10 1.185828e-74

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\theta_3$ (ma3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

$\phi_{2s}$ (sar2) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

Todos los parámetros son significativos.

Modelo 3: SARIMA(3,0,0)(2,1,2)

modelo3 <- Arima(data_T_DS, order = c(3,0,0), seasonal = c(2,0,2)) 
summary(modelo3)

Series: data_T_DS 
ARIMA(3,0,0)(2,0,2)[12] with non-zero mean 

Coefficients:
         ar1     ar2     ar3    sar1     sar2     sma1    sma2    mean
      0.0126  0.2686  0.3076  0.3996  -0.3256  -1.0836  0.3215  0.1507
s.e.  0.0735  0.0679  0.0768  0.4271   0.1387   0.4345  0.4177  0.0040

sigma^2 = 0.006116:  log likelihood = 214.16
AIC=-410.33   AICc=-409.34   BIC=-381.01

Training set error measures:
                      ME       RMSE        MAE       MPE     MAPE      MASE
Training set 0.001976132 0.07656065 0.05807014 -26.04676 77.44249 0.5075737
                   ACF1
Training set 0.02489652

matriz <- cbind(modelo3$coef, sqrt(diag(modelo3$var.coef)))
Z = matriz[,1]/matriz[,2]
valor_p <- 2*pnorm(-abs(Z))
valor_p

         ar1          ar2          ar3         sar1         sar2         sma1 
8.637860e-01 7.579762e-05 6.257770e-05 3.494727e-01 1.892622e-02 1.262932e-02 
        sma2    intercept 
4.414670e-01 0.000000e+00

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

$\phi_{2s}$ (sar2) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\theta_{2s}$ (sma2) ⟶ Como el valor p es mayor que la significancia $\alpha = 0.05$, se concluye que no es significativo.

De este modelo se concluye que $\theta_{2s}$ no es significativo. El modelo que resulta tras la eliminación de este coeficiente es el mismo hallado por el auto.arima, por lo tanto no se propone como nuevo modelo

Modelo 4: SARIMA(3,0,0)(2,1,0)

modelo4 <- Arima(data_T, order = c(3,0,0), seasonal = c(2,0,0))
summary(modelo4)

Series: data_T 
ARIMA(3,0,0)(2,0,0)[12] with non-zero mean 

Coefficients:
         ar1     ar2     ar3    sar1   sar2    mean
      0.1599  0.4136  0.3822  0.6286  0.318  2.5972
s.e.  0.0703  0.0599  0.0706  0.0786  0.082  1.3331

sigma^2 = 0.008927:  log likelihood = 179.87
AIC=-345.74   AICc=-345.16   BIC=-322.51

Training set error measures:
                     ME      RMSE        MAE       MPE     MAPE      MASE
Training set 0.01040882 0.0930814 0.07194017 0.2992048 2.916265 0.4356852
                    ACF1
Training set -0.01994847

matriz <- cbind(modelo4$coef, sqrt(diag(modelo4$var.coef)))
Z = matriz[,1]/matriz[,2]
valor_p <- 2*pnorm(-abs(Z))
valor_p

         ar1          ar2          ar3         sar1         sar2    intercept 
2.303223e-02 5.075962e-12 6.070740e-08 1.283519e-15 1.062447e-04 5.138602e-02

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\phi_{2s}$ (sar2) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

De este modelo se concluye que todos los parámetros son significativos.

Modelo 5: SARIMA(3,0,0)(2,1,1)

modelo5 <- auto.arima(data_T, approximation = FALSE)
summary(modelo5)

Series: data_T 
ARIMA(3,0,0)(2,1,1)[12] with drift 

Coefficients:
         ar1     ar2     ar3    sar1     sar2     sma1   drift
      0.0008  0.2681  0.3164  0.0645  -0.2576  -0.7389  0.0126
s.e.  0.0720  0.0678  0.0754  0.1187   0.1019   0.1041  0.0003

sigma^2 = 0.006082:  log likelihood = 214.12
AIC=-412.23   AICc=-411.44   BIC=-386.17

Training set error measures:
                      ME       RMSE        MAE        MPE     MAPE      MASE
Training set 0.001769708 0.07426635 0.05501052 0.07217916 2.086989 0.3331556
                   ACF1
Training set 0.02815937

matriz <- cbind(modelo5$coef, sqrt(diag(modelo5$var.coef)))
Z = matriz[,1]/matriz[,2]
valor_p <- 2*pnorm(-abs(Z))
valor_p

         ar1          ar2          ar3         sar1         sar2         sma1 
9.916184e-01 7.697476e-05 2.722135e-05 5.866133e-01 1.148761e-02 1.244180e-12 
       drift 
0.000000e+00

Como estos modelos son jerárquicos, es suficiente que miremos los coeficientes de órden superior para cada parámetro. Si este es significativo, los otros no podrán ser eliminados del modelo.

$\phi_3$ (ar3) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

$\phi_{2s}$ (sar2) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.
$\theta_{1s}$ (sma1) ⟶ Como el valor p es menor que la significancia $\alpha = 0.05$, se concluye que es significativo.

De este modelo se concluye que todos los parámetros son significativos.

Comparación modelos

Ahora se procede a comparar los modelos cuyos coeficientes son significativos. Usaremos como medidas de desempeño el Mean Absolute Percentual Error (MAPE), el criterio de información de Akaike AIC y el criterio bayesiano de Schwarz BIC.

MAPEs = cbind(MAPEmod1_2 = mape(modelo1_2$x, modelo1_2$fitted),
              MAPEmod2 = mape(modelo2$x, modelo2$fitted),
              MAPEmod4 = mape(modelo4$x, modelo4$fitted),
              MAPEmod5 = mape(modelo5$x, modelo5$fitted))
MAPEs

     MAPEmod1_2  MAPEmod2   MAPEmod4   MAPEmod5
[1,]  0.9001941 0.8210343 0.02916265 0.02086989

AICs=cbind(AICmod1_2 = AIC(modelo1_2), AICmod2 = AIC(modelo2),
           AICmod4 = AIC(modelo4), AICmod5 = AIC(modelo5))
AICs

     AICmod1_2   AICmod2  AICmod4   AICmod5
[1,]  -366.766 -396.3251 -345.736 -412.2311

BICs=cbind(BICmod1_2 = BIC(modelo1_2), BICmod2 = BIC(modelo2),
           BICmod4 = BIC(modelo4), BICmod5 = BIC(modelo5))
BICs

     BICmod1_2   BICmod2   BICmod4   BICmod5
[1,] -337.4485 -363.7501 -322.5092 -386.1711

Por unanimidad el mejor modelo es el modelo 5, pues para todas las medidas de desempeño tenidas en cuenta es el que obtuvo un valor más bajo. Esto nos asegura que sus predicciones son más confiables.

Pronósticos

Pronosticamos con un horizonte de 12 meses y un nivel de confianza del 95%. Nuestro modelo parece captar bien la esencia de la serie transformada y los pronósticos parecen acertados.

autoplot(forecast(modelo5, h=12, level=0.95))

Validación de supuestos

Para las pruebas estadísticas se tendrá en cuenta un nivel de significancia $\alpha=0.01$

checkresiduals(modelo5, test = FALSE)

Los residuales parecen ruido blanco.

Supuesto de normalidad

gghistogram(modelo5$residuals)
boxplot(modelo5$residuals)
qqPlot(modelo5$residuals)

[1] 185 200

Las pruebas gráficas nos muestran que, a pesar de que en principio parece haber una normalidad, en el qqplot se sugiere un problema de colas pesadas, pues estas se salen de los valores teóricos de normalidad.

\[ H_0:Los\;residuales\;siguen\;una\;distribución\;normal \]

\[ H_1:Los\;residuales\;no\;siguen\;una\;distribución\;normal \]

shapiro.test(modelo5$residuals)


    Shapiro-Wilk normality test

data:  modelo5$residuals
W = 0.98292, p-value = 0.01413

ad.test(modelo5$residuals)


    Anderson-Darling normality test

data:  modelo5$residuals
A = 1.1022, p-value = 0.006742

cvm.test(modelo5$residuals)


    Cramer-von Mises normality test

data:  modelo5$residuals
W = 0.19377, p-value = 0.006371

Los residuales del modelo pasan únicamente la prueba de normalidad de Shapiro-Wilk.

Supuesto de homocedasticidad
```
ggAcf(modelo5$residuals)
ggPacf(modelo5$residuals)
```
Gráficamente los residuales parecen tener varianza constante, solo una línea alrededor del rezago 18 se sale de las bandas para el caso de la ACF.

\[ H_0:La\;varianza\;de\;los\;residuales\;es\;constante\; \]

\[ H_1:La\;varianza\;de\;los\;residuales\;no\;es\;constante\; \]
```
bptest(modelo5$fitted~modelo5$residuals)
```
```
    studentized Breusch-Pagan test

data:  modelo5$fitted ~ modelo5$residuals
BP = 4.7638, df = 1, p-value = 0.02906
```
El valor p es mayor que el nivel de significancia, por lo tanto, no se rechaza la hipótesis nula. Se acepta el supuesto de homocedasticidad y se concluye que la varianza de los residuales es constante.
Supuesto de independencia

autoplot(modelo5$residuals)
ggAcf(modelo5$residuals)

La gráfica muestra una aparente independencia entre los residuales.

\[ H_0:Los\;residuales\;son\;independientes\; \]

\[ H_1:Los\;residuales\;no\;son\;independientes \]

Box.test(modelo5$residuals)


    Box-Pierce test

data:  modelo5$residuals
X-squared = 0.16176, df = 1, p-value = 0.6875

El valor p es mayor que el nivel de significancia. Se acepta la hipótesis nula y se valida el supuesto de independencia de los residuales.

Ten en cuenta:

A pesar de que en este caso no se hizo validación cruzada, esta estrategia es muy recomendada en el ajuste de modelos.