Aplicacion

Contexto de los datos

Los datos se encuentran en chickwts de R-project.Se realizó un experimento para medir y comparar la efectividad de varios,complementos alimenticios en la tasa de crecimiento de los pollos.Los pollitos recién nacidos se distribuyeron aleatoriamente en seis grupos,y cada grupo recibió un suplemento alimenticio diferente.Sus pesos en gramos después de seis semanas se dan junto con los tipos de alimento.El marco de datos contiene 71 observaciones sobre las siguientes 2 variables.

weight

Es una variable numérica que da el peso del pollo(gramos)

feed

Es un factor que da el tipo de alimentación.(horsebean,linseed,soybean,sunflower)

Seleccion del modelo

set.seed(31)
data(chickwts)
ts.plot(chickwts$weight, ylab="Peso en gramos", main="Evolucion del peso");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 tampoco parece ser 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.

w=log(chickwts$weight)
d<-diff(w)
ts.plot(d,ylab = "",main= "serie diferenciada");grid(nx = NULL, ny = NULL,lty = 2,col = "gray",lwd = 2)

# ACF y PACF muestrales de la serie original

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

#La ACF, presenta un decaimiento exponencial,por otra parte la PACF, se corta 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,se corta despues del rezago 1.La PACF presenta un decaimiento exponencial.

# Identificacion usando criterios de informacion
# entre los ordenes maximos de los polinomios

p=1
q=1
P=1
Q=1

# Entre los ordenes de diferenciacion

d=1
D=1
(maxfilas=(p+1)*(q+1)*(P+1)*(Q+1))

## [1] 16

ic_mod=matrix(rep(-99, times=maxfilas*6), nrow=maxfilas, ncol=6)
k=1
for(i in 0:p) {
  for(j in 0:q) {
    for(s in 0:P) {
      for(m in 0:Q) {
        modelo=Arima(w, order = c(i, d, j),
                     seasonal = list(order = c(s, D, m)), method = c("CSS-ML"))
        
        numpar=length(coef(modelo))
        T=length(residuals(modelo))
        ic_mod[k,1]=i
        ic_mod[k,2]=j
        ic_mod[k,3]=s
        ic_mod[k,4]=m
        ic_mod[k,5]=-2*(modelo$loglik/T)+(2*numpar)/T      # AIC
        ic_mod[k,6]=-2*(modelo$loglik/T)+(numpar*log(T))/T  # SIC
        k=k+1
      }
    }
  }
}
p_=ic_mod[,1]
q_=ic_mod[,2]
P_=ic_mod[,3]
Q_=ic_mod[,4]
AIC=ic_mod[,5]
SIC=ic_mod[,6]

# Tabla con ordenes de los modelos ajustados y los criterios de informacion
# ---------------------------------------------------------------------------

(Crit_Inf=cbind(p_, q_, P_, Q_, AIC, SIC))

##       p_ q_ P_ Q_       AIC       SIC
##  [1,]  0  0  0  0 1.2960130 1.2960130
##  [2,]  0  0  0  1 0.4739992 0.5058679
##  [3,]  0  0  1  0 0.9560883 0.9879570
##  [4,]  0  0  1  1 0.4309999 0.4947373
##  [5,]  0  1  0  0 0.4739992 0.5058679
##  [6,]  0  1  0  1 0.3305759 0.3943134
##  [7,]  0  1  1  0 0.4309999 0.4947373
##  [8,]  0  1  1  1 0.2638631 0.3594693
##  [9,]  1  0  0  0 0.9560883 0.9879570
## [10,]  1  0  0  1 0.4309999 0.4947373
## [11,]  1  0  1  0 0.9200849 0.9838224
## [12,]  1  0  1  1 0.4515260 0.5471322
## [13,]  1  1  0  0 0.4309999 0.4947373
## [14,]  1  1  0  1 0.2638631 0.3594693
## [15,]  1  1  1  0 0.4515260 0.5471322
## [16,]  1  1  1  1 0.2883688 0.4158437

#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
## 14  1  1  0  1 0.2638631 0.3594693
## 8   0  1  1  1 0.2638631 0.3594693
## 16  1  1  1  1 0.2883688 0.4158437
## 6   0  1  0  1 0.3305759 0.3943134
## 7   0  1  1  0 0.4309999 0.4947373
## 4   0  0  1  1 0.4309999 0.4947373
## 10  1  0  0  1 0.4309999 0.4947373
## 13  1  1  0  0 0.4309999 0.4947373
## 15  1  1  1  0 0.4515260 0.5471322
## 12  1  0  1  1 0.4515260 0.5471322
## 2   0  0  0  1 0.4739992 0.5058679
## 5   0  1  0  0 0.4739992 0.5058679
## 11  1  0  1  0 0.9200849 0.9838224
## 3   0  0  1  0 0.9560883 0.9879570
## 9   1  0  0  0 0.9560883 0.9879570
## 1   0  0  0  0 1.2960130 1.2960130

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

##    p_ q_ P_ Q_       AIC       SIC
## 14  1  1  0  1 0.2638631 0.3594693
## 8   0  1  1  1 0.2638631 0.3594693
## 6   0  1  0  1 0.3305759 0.3943134
## 16  1  1  1  1 0.2883688 0.4158437
## 7   0  1  1  0 0.4309999 0.4947373
## 4   0  0  1  1 0.4309999 0.4947373
## 10  1  0  0  1 0.4309999 0.4947373
## 13  1  1  0  0 0.4309999 0.4947373
## 2   0  0  0  1 0.4739992 0.5058679
## 5   0  1  0  0 0.4739992 0.5058679
## 15  1  1  1  0 0.4515260 0.5471322
## 12  1  0  1  1 0.4515260 0.5471322
## 11  1  0  1  0 0.9200849 0.9838224
## 3   0  0  1  0 0.9560883 0.9879570
## 9   1  0  0  0 0.9560883 0.9879570
## 1   0  0  0  0 1.2960130 1.2960130

#De acuerdo a la tabla,podemos observar que el modelo 14 presenta el menor criterio de informacion, por tanto es el mas adecuado.Es decir un ARIMA(1,1,1)(0,1,1)

#Especificacion y estimacion del modelo

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

## 
## Call:
## arima(x = w, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12), 
##     method = c("CSS-ML"))
## 
## Coefficients:
##          ar1      ma1     sma1
##       0.2588  -0.9087  -0.9995
## s.e.  0.1425   0.0730   0.8839
## 
## sigma^2 estimated as 0.0624:  log likelihood = -13.32,  aic = 34.64

Diagnostico residual

#Chequeo si los residuales son ruido blanco y posibles observaciones extremas

tsdiag(mod1)

#Otra forma de verificar si los residuales son ruido blanco

resmod1=residuals(mod1)
par(mfrow=c(2,1)) 
acf(resmod1)
pacf(resmod1)

#Los residuales no estan autocorrelacionados.Es decir son ruido blanco.

# Prueba de Ljung-Box

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

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

#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

shapiro.test(res_est)                #Prueba de Shapiro-Wilks

## 
##  Shapiro-Wilk normality test
## 
## data:  res_est
## W = 0.97816, p-value = 0.2516

# 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")

# Pronosticos para la serie transformada

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

## $pred
## Time Series:
## Start = 72 
## End = 95 
## Frequency = 1 
##  [1] 5.903647 5.985758 5.891002 5.955214 6.009289 5.898335 5.731113 5.873174
##  [9] 5.787306 5.751812 5.860209 5.925567 6.026918 6.109764 6.015198 6.079460
## [17] 6.133547 6.022596 5.855376 5.997437 5.911568 5.876075 5.984472 6.049830
## 
## $se
## Time Series:
## Start = 72 
## End = 95 
## Frequency = 1 
##  [1] 0.2739206 0.2886902 0.2929323 0.2953294 0.2973204 0.2992035 0.3010508
##  [8] 0.3028811 0.3047002 0.3065129 0.3083335 0.3102094 0.3207348 0.3239919
## [15] 0.3264761 0.3287184 0.3308895 0.3330321 0.3351574 0.3372689 0.3393688
## [22] 0.3414624 0.3435688 0.3457555

# Grafico
#plot(w,xlim=c(1949,1962), ylim=c(4,7))
#w.pred<-predict(mod1,n.ahead=24)
#lines(w.pred$pred,col="blue")
#lines(w.pred$pred+2*pasajt.pred$se,col="red",lty=3)
#lines(w.pred$pred-2*pasajt.pred$se,col="red",lty=3)

# Intervalos de prediccion para la serie transformada

pronost=w.pred$pred
pronost

## Time Series:
## Start = 72 
## End = 95 
## Frequency = 1 
##  [1] 5.903647 5.985758 5.891002 5.955214 6.009289 5.898335 5.731113 5.873174
##  [9] 5.787306 5.751812 5.860209 5.925567 6.026918 6.109764 6.015198 6.079460
## [17] 6.133547 6.022596 5.855376 5.997437 5.911568 5.876075 5.984472 6.049830

stdev=w.pred$se
stdev

## Time Series:
## Start = 72 
## End = 95 
## Frequency = 1 
##  [1] 0.2739206 0.2886902 0.2929323 0.2953294 0.2973204 0.2992035 0.3010508
##  [8] 0.3028811 0.3047002 0.3065129 0.3083335 0.3102094 0.3207348 0.3239919
## [15] 0.3264761 0.3287184 0.3308895 0.3330321 0.3351574 0.3372689 0.3393688
## [22] 0.3414624 0.3435688 0.3457555

liminft=pronost-2*stdev
limsupt=pronost+2*stdev

# Pronosticos e intervalos de prediccion para la serie original
# con correccion por sesgos de retransformacion

pronos_orig=exp(pronost+.5*stdev^2)
pronos_orig

## Time Series:
## Start = 72 
## End = 95 
## Frequency = 1 
##  [1] 380.3771 414.6475 377.6269 402.9546 425.5950 381.1131 322.6052 372.0554
##  [9] 341.6295 329.8993 367.8752 392.9500 436.3104 474.4942 432.0285 461.0413
## [17] 487.0129 436.1784 369.2748 425.9458 391.1745 377.8032 421.3616 450.1602

liminf=exp(liminft+.5*stdev^2)
limsup=exp(limsupt+.5*stdev^2)

#ts.plot(chickwts$weight, pronos_orig, liminf, limsup) 
#lines(chickwts$weight,col="black")
#lines(pronos_orig,col="blue")
#lines(liminf,col="red",lty=3)
#lines(limsup,col="red",lty=3)

# Listado de los pron?sticos e intervalos de prediccion

pronosticos=cbind(pronos_orig, liminf, limsup)
pronosticos

## Time Series:
## Start = 72 
## End = 95 
## Frequency = 1 
##    pronos_orig   liminf   limsup
## 72    380.3771 219.9328 657.8680
## 73    414.6475 232.7695 738.6389
## 74    377.6269 210.1963 678.4233
## 75    402.9546 223.2216 727.4045
## 76    425.5950 234.8266 771.3397
## 77    381.1131 209.4928 693.3279
## 78    322.6052 176.6778 589.0616
## 79    372.0554 203.0151 681.8468
## 80    341.6295 185.7360 628.3687
## 81    329.8993 178.7095 608.9970
## 82    367.8752 198.5571 681.5780
## 83    392.9500 211.2967 730.7717
## 84    436.3104 229.7253 828.6712
## 85    474.4942 248.2076 907.0823
## 86    432.0285 224.8738 830.0151
## 87    461.0413 238.9014 889.7357
## 88    487.0129 251.2659 943.9466
## 89    436.1784 224.0764 849.0479
## 90    369.2748 188.9016 721.8778
## 91    425.9458 216.9733 836.1849
## 92    391.1745 198.4260 771.1565
## 93    377.8032 190.8425 747.9215
## 94    421.3616 211.9507 837.6739
## 95    450.1602 225.4487 898.8484

MAPE

#Calcular MAPE

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] "11.0477382881477 %"

Conclusion