Aplicacion

Contexto de los datos

Aves de corral (pollo), precio al contado de aves enteras, muelles de Georgia, centavos de dólar EE.UU. por libra

Grafica del precio del pollo en centavos por libra entre los años 2005 y 2015

ts.plot(chicken, ylab="Centavos por libra", main="Precios del pollo");grid(nx = NULL, ny = NULL,lty = 2,col = "gray",lwd = 2)

La serie no es estacionaria.Se puede apreciar que no esta en torno una media,y la varianza parece incrementar con el tiempo,es decir no estacionaria en varianza.Al aplicar una transformacion logaritmica se puede corregir el problema de estacionariedad en varianza,por otro lado para corregir la no estacionariedad en media procedemos a diferenciar la serie.Para tener una serie estacionaria en media y varianza.

Prueba de Dickey-Fuller

La Prueba de Dickey-Fuller busca determinar la existencia o no de raices unitarias en una serie de tiempo La hipótesis nula de esta prueba es que existe una raíz unitaria en la serie.

Hipotesis:

Ho:La serie no es estacionaria

Ha:La series es estacionaria

adf.test(chicken)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  chicken
## Dickey-Fuller = -2.3955, Lag order = 5, p-value = 0.4109
## alternative hypothesis: stationary

A un nivel de significancia de 0.05 no se rechaza Ho, es decir la serie no es estacionaria.

Grafica de la serie diferenciada

ts.plot(d,ylab = "",main= "Serie diferenciada");grid(nx = NULL, ny = NULL,lty = 2,col = "gray",lwd = 2)

Ahora se puede apreciar que esta en torno a una media no depende del tiempo, y la varianza tampoco incrementa en el tiempo.

Graficos serie original y en logaritmos

par(mfrow=c(2,1))
plot(chicken, main="Serie Original") 
plot(log(chicken), main="Logaritmo de la serie")

Comparacion entre la serie y al aplicar una transformacion logaritimica.

Descomposicion de la serie

autoplot(stl(chicken, s.window = "periodic"), ts.colour = "blue")

Tendencia: Tiene una tendencia creciente, es decir el precio aumenta con el paso del tiempo

Estacionalidad:Se presenta un comportamiento de estacionalidad cada año.

ACF y PACF muestrales de la serie original

par(mfrow=c(2,1)) 
acf(chicken, lag.max=48)     # ACF hasta el rezago 48
pacf(chicken, lag.max=48)    # PACF hasta el rezago 48

La ACF presenta un decaimiento exponencial, y la PACF un corte despues del rezago 1.

ACF y PACF muestrales de la primera diferencia ordnaria(1-B).

par(mfrow=c(2,1)) 
acf(diff(w), lag.max=48)     # ACF hasta el rezago 48
pacf(diff(w), lag.max=48)    # PACF hasta el rezago 48

La ACF, presenta un decaimiento exponencial como ondas amortiguadas y la PACF un corte despues del rezago 2, es decir un AR(2) y un MA(0),Lo que podria sugerir un ARIMA(2,1,0).

ACF y PACF muestrales estacionales

par(mfrow=c(2,1)) 
acf(w, lag=12, lag.max=48)     # ACF hasta el rezago 48
pacf(w, lag=12, lag.max=48)    # PACF hasta el rezago 48

#acf(diff(w), lag=12, lag.max=48)     # ACF hasta el rezago 48
#pacf(diff(w), lag=12, lag.max=48)    # PACF hasta el rezago 48

La ACF, presenta un decaimiento exponencial, y su PACF, se corta despues del rezago 1,es decir un AR(1) y MA(0) lo que podria sugerir un ARMA(1,0) o ARIMA(1,0,0)

Modelo

Ordenamiento de los criterios de informacion e menor a mayor

Crit_Inf=data.frame(Crit_Inf)
(aic=Crit_Inf[order(Crit_Inf$AIC),])

##   p_ q_ P_ Q_       AIC       SIC
## 6  2  0  1  0 -6.811207 -6.757991
## 4  1  0  1  0 -6.751626 -6.716148
## 5  2  0  0  0 -6.701521 -6.666043
## 3  1  0  0  0 -6.603194 -6.585455
## 2  0  0  1  0 -5.972472 -5.954733
## 1  0  0  0  0 -5.834296 -5.834296

(sic=Crit_Inf[order(Crit_Inf$SIC),])

##   p_ q_ P_ Q_       AIC       SIC
## 6  2  0  1  0 -6.811207 -6.757991
## 4  1  0  1  0 -6.751626 -6.716148
## 5  2  0  0  0 -6.701521 -6.666043
## 3  1  0  0  0 -6.603194 -6.585455
## 2  0  0  1  0 -5.972472 -5.954733
## 1  0  0  0  0 -5.834296 -5.834296

De acuerdo a la tabla,podemos observar que el modelo que presenta el menor criterio de informacion, por tanto es el mas adecuado.Es decir un SARIMA(2,1,0)*(1,0,0).SARIMA combina una parte estacional y no estacional con S=12

Especificacion y estimacion del modelo

(mod1=arima(w, order = c(2, 1, 0), seasonal = list(order = c(1, 0, 0), period = 12),method = c("CSS-ML")))

## 
## Call:
## arima(x = w, order = c(2, 1, 0), seasonal = list(order = c(1, 0, 0), period = 12), 
##     method = c("CSS-ML"))
## 
## Coefficients:
##          ar1      ar2    sar1
##       0.9419  -0.2682  0.3469
## s.e.  0.0733   0.0738  0.0712
## 
## sigma^2 estimated as 5.92e-05:  log likelihood = 616.01,  aic = -1226.02

Diagnostico residual

tsdiag(mod1)

Los residuales no estan autocorrelacionados.Es decir son ruido blanco,se encuentran dentro de la franja.Tambien son estadisticamente diferentes de cero.

Prueba de Ljung-Box

Box.test(resmod1, lag = 36, type = "Ljung-Box")  #¿Porque lag=36?

## 
##  Box-Ljung test
## 
## data:  resmod1
## X-squared = 41.85, df = 36, p-value = 0.2318

Los residuales se distribuyen de forma independiente a un nivel de significancia de 0.05

Chequeo de normalidad(Grafico cuantil-cuantil)

res_est=mod1$residuals
qqnorm(res_est,  xlab = "Cuantiles Teoricos", ylab = "Cuantiles Muestrales")
qqline(res_est)

Podemos observar un buen ajuste a la recta.Los residuales parecen seguir una distribucion normal.

Prueba de normalidad(Prueba de Shapiro-Wilks)

shapiro.test(res_est)

## 
##  Shapiro-Wilk normality test
## 
## data:  res_est
## W = 0.99277, p-value = 0.5143

No se puede rechazar normalidad a un nivel de significancia de 0.05.Es decir los residuales siguen una distribucion normal.

Pronostico

Valores ajustados

ajust=w-resmod1

Grafico para los valores ajustados

ts.plot(w, ajust)
lines(w, col="black")
lines(ajust, col="red")

Se puede observar muy buen ajuste, la diferencia entre lo observado y lo predicho es “minima”

Pronosticos para la serie transformada

(w.pred=predict(mod1, n.ahead = 24, se.fit = TRUE))#¿Por que 24?

## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2016                                                                4.709790
## 2017 4.698042 4.696768 4.695116 4.695101 4.696443 4.696013 4.694834 4.693490
## 2018 4.689414 4.688972 4.688399 4.688394 4.688860 4.688710 4.688301         
##           Sep      Oct      Nov      Dec
## 2016 4.706830 4.703823 4.700720 4.699705
## 2017 4.692463 4.691419 4.690343 4.689991
## 2018                                    
## 
## $se
##              Jan         Feb         Mar         Apr         May         Jun
## 2016                                                                        
## 2017 0.047695892 0.053297544 0.058342554 0.062961057 0.067247302 0.071267756
## 2018 0.105227312 0.109936790 0.114447578 0.118781356 0.122958436 0.126995974
##              Jul         Aug         Sep         Oct         Nov         Dec
## 2016             0.007694249 0.016806183 0.025897950 0.034141260 0.041370604
## 2017 0.075069418 0.079526868 0.084568686 0.089872030 0.095165956 0.100302444
## 2018 0.130908030

Grafico con horizonte

ts.plot(chicken, pronos_orig, liminf, limsup) 
lines(chicken,col="black")
lines(pronos_orig,col="blue")
lines(liminf,col="red",lty=3)
lines(limsup,col="red",lty=3)

Listado de los pronosticos e intervalos de prediccion

pronosticos=cbind(pronos_orig, liminf, limsup)
head(pronosticos)

##          pronos_orig    liminf   limsup
## Aug 2016    111.0321 109.33660 112.7540
## Sep 2016    110.7163 107.05673 114.5010
## Oct 2016    110.4053 104.83234 116.2745
## Nov 2016    110.0905 102.82417 117.8704
## Dec 2016    110.0088 101.27295 119.4983
## Jan 2017    109.8570  99.86183 120.8525

MAPE(Error Porcentual Absoluto Medio)

split = sample.split(w, SplitRatio = 0.9)
train = subset(w, split == TRUE)
test1 = subset(w, split == FALSE)
mp= as.character(MAPE(pronost,test1)*100)
paste(mp,"%")

## [1] "7.2402314736348 %"

Conclusion

De acuerdo a los criterio de informacion un (ARIMA(2,1,0)(1,0,0)),tiene el menor criterio de informacion.Sus diagnosticos residuales son buenos, es decir cumple con los supuestos del modelo.Ademas el MAPE es 7.2402(Bueno).Por lo tanto el modelo escogido tiene un buen desempeño. Es decir menor criterio de informacion,cumple los supuestos y tiene buen pronostico.

Aplicacion

Daniel Vargas Zapata