Quiz 3

Authors

Kevin Duvan Barrera Pezca

Mateo Hasane Vega García

Se va a hacer un análisis de outliers para la serie de pasajeros, primero se hace un análisis descriptivos de la serie.

data(AirPassengers)
plot(AirPassengers)

Con esta gráfica se sospecha que la serie tiene estacionalidad y se observa que tiene varianza no constante y tendencia, se hace una transformación logarítimica para estabilizar la varianza.

lAirpassengers=log(AirPassengers)
plot(lAirpassengers)

Ahora se hace la prube a de Dickey - Fuller para detectar estacionalidad

tseries::adf.test(lAirpassengers,k=12)


    Augmented Dickey-Fuller Test

data:  lAirpassengers
Dickey-Fuller = -1.5325, Lag order = 12, p-value = 0.7711
alternative hypothesis: stationary

Como el p valor es mayor al nivel de significancia no se rechaza la hipótesis nula, por lo tanto hay suficiente evidencia de que los datos son no estacionarios, se ajusta un modelo SARIMA.

dlAirpassengers=diff(lAirpassengers,lag=1)
par(mfrow=c(2,1))
plot(lAirpassengers, ylab = 'datos transformados')
plot(dlAirpassengers, ylab = 'datos diferenciados')

tseries::adf.test(dlAirpassengers,k = 15)


    Augmented Dickey-Fuller Test

data:  dlAirpassengers
Dickey-Fuller = -3.3376, Lag order = 15, p-value = 0.0678
alternative hypothesis: stationary

Al hacer el test nuevamente no se rechaza la hipótesis nula, todavía se debería hacer otra diferenciación.

monthplot(dlAirpassengers)

Con este gráfico se observa que todavía existe estacionalidad, se ejecuta nsdiff() para revisar si los datos requieren otra diferenciación todavía,

nsdiffs(dlAirpassengers)

[1] 1

El resultado muestra que todavía falta una diferenciación, una vez estabilizada la varianza y eliminada la estacionariedad, se hace un gráfico de autocorrelación y autocorrelación parcial.

DdlAirpassengers=diff(dlAirpassengers,lag=12)
acf(dlAirpassengers)

pacf(dlAirpassengers,lag.max = 48)

Ajuste del Modelo

Dados los gráficos anteriores se escoge un modelo SARIMA(2, 1, 3)(0, 1, 1)

modelo = Arima(AirPassengers, c(1, 1, 4),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0) 
# 1 diferenciación y transformación logarítmica
coeftest(modelo)


z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1   0.944009   0.048220 19.5773 < 2.2e-16 ***
ma1  -1.373330   0.147916 -9.2845 < 2.2e-16 ***
ma2   0.447062   0.165070  2.7083  0.006763 ** 
ma3  -0.222342   0.184917 -1.2024  0.229214    
ma4   0.148623   0.104568  1.4213  0.155228    
sma1 -0.553486   0.077168 -7.1725 7.365e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

modelo

Series: AirPassengers 
ARIMA(1,1,4)(0,1,1)[12] 
Box Cox transformation: lambda= 0 

Coefficients:
         ar1      ma1     ma2      ma3     ma4     sma1
      0.9440  -1.3733  0.4471  -0.2223  0.1486  -0.5535
s.e.  0.0482   0.1479  0.1651   0.1849  0.1046   0.0772

sigma^2 = 0.00134:  log likelihood = 247.06
AIC=-480.11   AICc=-479.2   BIC=-459.98

Se quitan los parámetros no significativos,

modeloref = Arima(AirPassengers, c(1, 1, 4),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0, fixed = c(NA, NA, NA, 0, 0, NA)) 
# 1 diferenciación y transformación logarítmica
coeftest(modeloref)


z test of coefficients:

      Estimate Std. Error  z value  Pr(>|z|)    
ar1   0.900474   0.090944   9.9014 < 2.2e-16 ***
ma1  -1.356290   0.105940 -12.8024 < 2.2e-16 ***
ma2   0.341780   0.111878   3.0549  0.002251 ** 
sma1 -0.552876   0.077086  -7.1722  7.38e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

modeloref

Series: AirPassengers 
ARIMA(1,1,4)(0,1,1)[12] 
Box Cox transformation: lambda= 0 

Coefficients:
         ar1      ma1     ma2  ma3  ma4     sma1
      0.9005  -1.3563  0.3418    0    0  -0.5529
s.e.  0.0909   0.1059  0.1119    0    0   0.0771

sigma^2 = 0.001295:  log likelihood = 246.02
AIC=-482.04   AICc=-481.56   BIC=-467.67

Análisis de residuales

res = modeloref$residuals
plot(res)

Los residuos no presentan tendencia y están alrededor de 0, sin embargo se observa que hay datos alejados de la media, estos podrían ser outliers.

acf(res)

acf(res^2)

Se observa que los residuos se comportan como un ruido blanco, es decir, se cumple nuestro supuesto.

Los residuos al cuadrado se comportan bien, es decir, no hay varianza condicional no constante.

Box.test(res, lag=12 , type = "Ljung-Box", fitdf = 2)


    Box-Ljung test

data:  res
X-squared = 8.7911, df = 10, p-value = 0.552

No se rechaza la hipótesis de no autocorrelación.

tseries::jarque.bera.test(res)


    Jarque Bera Test

data:  res
X-squared = 11.501, df = 2, p-value = 0.003182

Se rechaza la hipótesis nula, no hay evidencia suficiente para asegurar que los errores provengan de una dsitribución normal.

hist(res)

plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")
lines(LI,type="S",col="red")

plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")                      
lines(LQS,type="S",col="red")                                                                           
lines(LQI,type="S",col="red")

Se observa que el gráfico cumsum al cuadrado se sale de los intervalos.

Outliers

coef= coefs2poly(modelo)
coef

$arcoefs
 [1]  1.9440085 -0.9440085  0.0000000  0.0000000  0.0000000  0.0000000
 [7]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  1.0000000
[13] -1.9440085  0.9440085

$macoefs
 [1] -1.37333015  0.44706158 -0.22234178  0.14862298  0.00000000  0.00000000
 [7]  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 -0.55348580
[13]  0.76011874 -0.24744224  0.12306302 -0.08226071

attr(,"class")
[1] "ArimaPars"

outliers= locate.outliers(res, coef)
outliers

  type ind     coefhat     tstat
1   AO  29  0.08213950  4.231479
2   AO  38  0.06854757  3.531259
3   AO  62 -0.07744486 -3.989274
4   AO 135 -0.09934798 -4.510079
5   LS  24  0.07523411  3.623006
6   LS  54 -0.08606390 -4.144449

Se encuentra un outlier de tipo aditivo en la observación 29, 38, 62 y de tipo cambio de nivel en las observacion 24 y 54.

n=length(AirPassengers)
xreg = outliers.effects(outliers, n)

modelo1 = Arima(AirPassengers, c(1, 1, 2),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0, xreg = xreg)
coeftest(modelo1)


z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
ar1   -0.465964   0.350626 -1.3289 0.1838655    
ma1    0.145309   0.345804  0.4202 0.6743348    
ma2   -0.036408   0.187405 -0.1943 0.8459598    
sma1  -0.448477   0.086118 -5.2077 1.912e-07 ***
AO29   0.090749   0.020400  4.4486 8.644e-06 ***
AO38   0.067881   0.020034  3.3883 0.0007032 ***
AO62  -0.066747   0.019994 -3.3384 0.0008426 ***
AO135 -0.103743   0.022575 -4.5955 4.317e-06 ***
LS24   0.081188   0.024084  3.3710 0.0007490 ***
LS54  -0.103198   0.023905 -4.3170 1.582e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

modelo1

Series: AirPassengers 
Regression with ARIMA(1,1,2)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1     ma1      ma2     sma1    AO29    AO38     AO62    AO135
      -0.4660  0.1453  -0.0364  -0.4485  0.0907  0.0679  -0.0667  -0.1037
s.e.   0.3506  0.3458   0.1874   0.0861  0.0204  0.0200   0.0200   0.0226
        LS24     LS54
      0.0812  -0.1032
s.e.  0.0241   0.0239

sigma^2 = 0.0008681:  log likelihood = 279.83
AIC=-537.66   AICc=-535.44   BIC=-506.03

modelo1ref = Arima(AirPassengers, c(1, 1, 2),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0, xreg = xreg, fixed = c(0, 0, 0, NA, NA, NA, NA, NA, NA, NA))
coeftest(modelo1ref)


z test of coefficients:

       Estimate Std. Error z value  Pr(>|z|)    
sma1  -0.446800   0.098109 -4.5541 5.261e-06 ***
AO29   0.104292   0.018250  5.7145 1.101e-08 ***
AO38   0.059098   0.018196  3.2479 0.0011627 ** 
AO62  -0.069066   0.018166 -3.8020 0.0001436 ***
AO135 -0.105709   0.021122 -5.0047 5.594e-07 ***
LS24   0.080245   0.026085  3.0762 0.0020963 ** 
LS54  -0.089974   0.025458 -3.5342 0.0004091 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

modelo1ref

Series: AirPassengers 
Regression with ARIMA(1,1,2)(0,1,1)[12] errors 
Box Cox transformation: lambda= 0 

Coefficients:
      ar1  ma1  ma2     sma1    AO29    AO38     AO62    AO135    LS24     LS54
        0    0    0  -0.4468  0.1043  0.0591  -0.0691  -0.1057  0.0802  -0.0900
s.e.    0    0    0   0.0981  0.0183  0.0182   0.0182   0.0211  0.0261   0.0255

sigma^2 = 0.0009409:  log likelihood = 272.99
AIC=-529.99   AICc=-528.81   BIC=-506.99

resi_analisis= residuals(modelo1ref)
coef_analisis= coefs2poly(modelo1ref)
outliers_analisis= locate.outliers(resi_analisis,coef_analisis)
outliers_analisis

[1] type    ind     coefhat tstat  
<0 rows> (or 0-length row.names)

Ya no hay outliers en el modelo.

res = modelo1ref$residuals
plot(res)

Box.test(res, lag=12 , type = "Ljung-Box", fitdf = 2)


    Box-Ljung test

data:  res
X-squared = 25.96, df = 10, p-value = 0.003795

tseries::jarque.bera.test(res)


    Jarque Bera Test

data:  res
X-squared = 0.43247, df = 2, p-value = 0.8055

No se rechaza la hipótesis nula, por lo tanto aun evindecia para segurar que los datos provienen de una distribución normal.

hist(res)

plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")
lines(LI,type="S",col="red")

plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")                      
lines(LQS,type="S",col="red")                                                                           
lines(LQI,type="S",col="red")

Se observa que el gráfico cumsum al cuadrado se comporta mejor después de la eliminación de los outliers.