library(tseries)
library(forecast)
library(TSA)
library(ggfortify)
library(ggplot2)
require(knitr)
require(pander)

J. Crew es una compañía de moda conocida por sus ventas a jóvenes profesionales, llegando a tener como su gran Cliente a Michelle Obama. Se usará los registros correspondientes a los ingresos trimestrales desde el primer trimestre de 2003 hasta el último trimestre de 2012. Los registros están en miles de dolares. Base de datos encontrada en DASL, software de análisis y exploración de datos.

Quarter Year Period Revenue Time
1 2003 0 161.5 2003.00
2 2003 1 167.1 2003.25
3 2003 2 151.0 2003.50
4 2003 3 210.4 2003.75
1 2004 4 145.7 2004.00
2 2004 5 188.6 2004.25

No se tomarán los últimos 3 registros correspondientes a los ingresos de la marca J. Crew de los últimos 3 trimestres del año 2012 con el fin de pronosticarlos y luego proceder a comparar los pronosticos con los registros reales con el uso de la validación cruzada.

Descripción

ingresos<-ts(jcrew[1:37,4],start = 2003,frequency = 4)
autoplot(ingresos, ts.colour = "turquoise",ts.size = 1.2,main = "Informe de Ingresos Trimestrales de J. Crew")

Según el gráfico presentado, se puede observar una tendencia en los ingresos, lo que quiere decir que a medida del tiempo la empresa J. Crew va aumentando sus ventas y por ende sus ingresos.Adicionalmente, se puede observar en nuestra serie una fuerte componente estacional, presentando picos altos en el cuarto trimestre de cada año y picos bajos en el 1er trimestre. De acuerdo a la varianza de la serie, podemos concluir en base a una interpretación visual, que la serie no presenta heterocedasticidad significativa.

Pasos de ajuste de modelos propuestos por Box- Jen-kins

Se tomará un nivel de significancia de 0.03

Paso 1. Transformaciones necesarias para lograr la estacionariedad

Criterio 1. Aplicar una prueba de raiz unitaria para determinar la cantidad de diferencias ordinarias.

\[ H_o: Hay \ raíz \ unitaria. \ (Serie \ no\ estacionaria) \\ H_a: No \ hay \ raíz \ unitaria. \ (Serie \ estacionaria) \]

Test de Dicket Fuller aplicada a la 1era diferencia ordinaria.

dtasa1<-diff(ingresos,lag=1,differences=1)
dtasa2<-diff(ingresos,lag=1,differences=2)
adf.test(dtasa1) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dtasa1
## Dickey-Fuller = -2.6991, Lag order = 3, p-value = 0.3011
## alternative hypothesis: stationary

Con un p-valor de 0.3011 se concluye que la serie diferenciada ordinariamente una vez no es estacionaria.

Test de Dicket Fuller aplicada a la 2da diferencia ordinaria.

adf.test(dtasa2) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dtasa2
## Dickey-Fuller = -3.8908, Lag order = 3, p-value = 0.02539
## alternative hypothesis: stationary

Con un p-valor de 0.02539 se concluye que la serie diferenciada ordinariamente dos vez sí es estacionaria.

Finalmente, este criterio nos propone diferenciar ordinariamente 2 veces la serie

Criterio 2: Método indireto, explorar el rápido decaimiento de la Función de autocorrelación.

acf(dtasa1, lag.max=40,main="FAC de la serie diferenciada una vez ordinariamente")

Según este criterio es suficiente con diferenciar ordinariamente una vez, pues su función de Autocorrelación decae rapidamente hacia cero. Adicionalmente, al parecer, es necesario alguna diferenciación estacional, porque la autocorrelación cada trimestre decae lentamente hacia cero.

Criterio 3: Elegir las diferencias que minimicen la desviación estándar (Propuesto por Victor Guerrero).La desviación disminuye hasta llegar a la estacionariedad y comienza a aumentar.

Diferenciación Ordinaria

s00=sd(ingresos)
s01=sd(diff(ingresos,lag=1,differences=1))
s10=sd(diff(ingresos,lag=1,differences=2))
kable(data.frame(Sx=round(c(s00,s01,s10),2)))
Sx
107.51
43.20
73.53

Según este criterio es suficiente con diferenciar una vez de forma ordinaria la serie para conseguir estacionariedad

Diferenciación Estacional

s0=sd(dtasa1)
s1=sd(diff(dtasa1,lag=4,differences=1))
s11=sd(diff(ingresos,lag=1,differences=2))
kable(data.frame(Sx=round(c(s0,s1,s11),2)))
Sx
43.20
21.89
73.53

Según este criterio es necesario diferenciar una vez de forma estaccional para conseguir estacionariedad.

Finalmente, se aplica una prueba de raiz unitaria a la serie con una diferencia ordinaria y una estacional.

dtasa<-diff(dtasa1,lag = 4,differences = 1)

adf.test(dtasa)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dtasa
## Dickey-Fuller = -4.0684, Lag order = 3, p-value = 0.01936
## alternative hypothesis: stationary
autoplot(dtasa, ts.colour = "turquoise",ts.size = 1.2,main = "Serie Estacionaria") + theme(plot.title = element_text(hjust = 0.5,face="bold",size=15))

En conclusión, de acuerdo a los resultados y la gráfica expuesta se concluye que la serie con una diferencia ordinaria y una estacional ya es estacionaria.

Paso 2. Identificación del modelo. (Mínimo 2)

opcion = par(mfrow=c(2,1))
acf(dtasa,lag.max=25) 
pacf(dtasa,lag.max=25)

En base a las gráficas de la Función de Autocorrelación y Autocorrelación Parcial, se propone dos modelos:

  • El primero con un AR(2) en parte Estacional.

\[ (1-\phi_1 B^4 -\phi_2 B^8)(1-B^4)(1-B)=a_{t} \]

  • El segundo con un MA(1) en parte Estacional.

\[ (1-B^4)(1-B)=(1-\theta_1B^4)a_t \]

Paso 3. Estimación de parámetros

Ahora se procede a estimar los parámetros de los dos modelos y presentar los coeficientes resultantes junto con sus intervalos de confianza.

1er modelo

modelo1=arima(ingresos,order=c(0,1,0),seasonal=list(order=c(2,1,0),period=4))
kable((round(confint(modelo1),2)))
2.5 % 97.5 %
sar1 -0.92 -0.19
sar2 -0.55 0.20

Algunos parametros no son significativos, y debido a que este modelos sí cumple con el supuesto de Normalidad,como se verá más adelante, entonces se quitarán aquellos parametros cuyo intervalo de confianza pase por cero, es decir el parámetro 2.

modelo1<-arima(ingresos,order=c(0,1,0),seasonal=list(order=c(2,1,0),period=4),fixed = c(NA,0))
kable(data.frame(Coef=round(modelo1$coef,4)))
Coef
sar1 -0.4924
sar2 0.0000 
kable(round(confint(modelo1),2))
2.5 % 97.5 %
sar1 -0.82 -0.16
sar2 NA NA

\[ (1+0.4924 B^4)(1-B^4)(1-B)=a_{t} \]

2do modelo

modelo2<-arima(ingresos,order = c(0,1,0),seasonal = list(order=c(0,1,1),period=4))
kable(data.frame(Coef=round(modelo2$coef,4)))
Coef
sma1 -0.5448 
kable(round(confint(modelo2),2))
2.5 % 97.5 %
sma1 -0.86 -0.23

\[ (1-B^4)(1-B)=(1+ 0.5448B^4)a_t \]

Paso 4. Verificación de supuestos de los residuales

Media cero

1er modelo

t.test(modelo1$residuals)
## 
##  One Sample t-test
## 
## data:  modelo1$residuals
## t = 0.96083, df = 36, p-value = 0.343
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -3.147850  8.815677
## sample estimates:
## mean of x 
##  2.833913

Con un p-valor de 0.343 se acepta la hipotesis nula lo que quiere decir que los residuales del primer modelo no difieren significativamente de cero.

2do modelo

t.test(modelo2$residuals)
## 
##  One Sample t-test
## 
## data:  modelo2$residuals
## t = 1.1567, df = 36, p-value = 0.255
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.507216  9.163521
## sample estimates:
## mean of x 
##  3.328153

con un p-valor de 0.255 se acepta la hipotesis nula lo que quiere decir que los residuales del segundo modelo no difieren significativamente de cero.

Normalidad

1er modelo

shapiro.test(modelo1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo1$residuals
## W = 0.93886, p-value = 0.04244

Con un p-valor de 0.04244 se concluye que los residuales del modelo 1 distribuyen Normal.

2do modelo

shapiro.test(modelo2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo2$residuals
## W = 0.97233, p-value = 0.4752

Con un p-valor de 0.4752 se concluye que los residuales del modelo 2 distribuyen Normal.

Blanquitud

1er modelo

Box.test(residuals(modelo1),lag=20,type="Ljung" )
## 
##  Box-Ljung test
## 
## data:  residuals(modelo1)
## X-squared = 19.92, df = 20, p-value = 0.4629
tsdiag(modelo1,gof.lag=10)

La estadística Q-Ljung Box arrojó un p-valor de 0.4629, por lo que se puede concluir que las primeras 20 autocorrelaciones de los residuales son simultaneamente cero. Además, Según la grafica de la Función de Autocorrelación se observa que las autocorrelaciones no se sale de las bandas. Así mismo, en el gráfico de los p-valores todos sobrepasan la bandita que representa el nivel de significancia (0.05).

2do modelo

Box.test(residuals(modelo2),lag=20,type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(modelo2)
## X-squared = 19.746, df = 20, p-value = 0.4739
tsdiag(modelo2,gof.lag=10)

Con un p-valor de 0.4739 arrojado por la estadística Q-Ljung Box se concluye que las primeras 20 autocorrelaciones de los residuales son simultaneamente cero. Además, según la grafica de la Función de Autocorrelación se observa que los valores no se salen de las bandas. Adicionalmente, se observa que en el gráfico de los p-valores todos sobrepasan la bandita que representa el nivel de significancia (0.05), por lo que se puede concluir blanquitud en los residuales del 2do modelo.

Paso 5. Elección del mejor modelo.

A continuación se verá los criterios de información de los dos modelos.

kma=length(modelo1$coef)
kar=length(modelo2$coef)
n = length(dtasa1)
n
## [1] 36
# AIC Returned Value
f<-log(modelo1$sigma2) + (n+2*kma)/n
g<-log(modelo2$sigma2) + (n+2*kar)/n

# AICc
a<-log(modelo1$sigma2) + (n+kma)/(n-kma-2) 
b<-log(modelo2$sigma2) + (n+kar)/(n-kar-2)

# BIC
d<-log(modelo1$sigma2) + kma*log(n)/n 
e<-log(modelo2$sigma2) + kar*log(n)/n

tab<-data.frame(Modelo=c("Modelo 1","Modelo 2"),
                AIC=c(f,g),AICc=c(a,b),BIC=c(d,e))
pander(tab)
Modelo AIC AICc BIC
Modelo 1 7.028 7.105 6.116
Modelo 2 6.934 7 5.978

Los 3 estadísticas seleccionan el modelo 2, por que presenta criterios de información menores.

Paso 6. Pronosticos de los modelos

A continuación, se predecirá las ventas de los 3 trimestres siguientes con los dos modelos ajustados y se comparará las predicciones con las observaciones reales. Esto con el fin de decidir cuál de los dos tiene mejor poder predictivo.

pron.m1 = predict(modelo1, n.ahead=3) 
pron.m2<-predict(modelo2,n.ahead = 3)

1er modelo

kable(data.frame(Trim=c("T2","T3","T4"),Pronosticos=round(pron.m1$pred,2),Reales=jcrew[38:40,4]))
Trim Pronosticos Reales
T2 513.29 525.5
T3 546.67 555.8
T4 592.25 642.9

2do modelo

kable(data.frame(Trim=c("T2","T3","T4"),Pronosticos=round(pron.m2$pred,2),Reales=jcrew[38:40,4]))
Trim Pronosticos Reales
T2 516.37 525.5
T3 552.90 555.8
T4 600.49 642.9

A continuación se presentará las estadísticas correspondientes al Error Cuadratico Medio y Error Porcentual Absoluto para concluir qué modelo predice mejor las ventas futuras.

pred1=pron.m1
pred2=pron.m2

observed=jcrew[38:40,4]

cme1=sqrt(mean((pred1$pred-observed)^2))
cme2=sqrt(mean((pred2$pred-observed)^2))

# MAPE es la variación porcentual
mape1=mean(abs((observed-pred1$pred)*100/observed))
mape2=mean(abs((observed-pred2$pred)*100/observed))

kable(data.frame(Modelo=c("uno","dos"),CME=round(c(cme1,cme2),2),MAPE=round(c(mape1,mape2),2)))
Modelo CME MAPE
uno 30.54 3.95
dos 25.10 2.95

El Error Cuadratico Medio mide el promedio de la distancia de lo predicho a lo observado, y el Error Porcentual Absoluto mide el promedio de la variación porcentual de lo predicho entre lo observado. En este orden de ideas, de acuerdo a los resultados de los criterios MAPE y ECM y teniendo como referencia las diferencias entre lo pronosticado por los dos modelos y los registros reales se llega a la conclusión que el mejor modelo en términos de poder predictivo es el modelo 2.

\[ (1-B^4)(1-B)=(1+ 0.5448B^4)a_t \]

A continuación se presenta la gráfica mostrando los pronosticos del modelo elegido junto con su intervalo de confianza, además se grafica los registros reales.

U = pron.m2$pred + 2*pron.m2$se
L = pron.m2$pred - 2*pron.m2$se

ingresos_real<-ts(jcrew[38:40,4],start = 2012.25,frequency = 4)

ts.plot(ingresos,pron.m2$pred, col=1:2, type="o",lwd=3)
lines(U, col="blue", lty="dashed",lwd=2)
lines(L, col="blue", lty="dashed",lwd=2)
lines(ingresos_real,col="green",lty=,lwd=2)
abline(v=2012.25,lty="dotted")
legend("topleft",legend = c("Registros Reales","Registros Pronosticados"),col=c("green","red"),lwd=c(2,2),cex=1.2, box.col="purple",box.lty=2, box.lwd=2)

Resumen

Inicialmente, de acuerdo al uso de los 3 criterios para hallar las diferencias necesarios para lograr la estacionariedad en la serie, se optó por diferenciar una vez tanto ordinaria como estacionalmente. Fue complicado proponer modelos de acuerdo a la grafica de la función de autocorrelación y la función de autocorrelación parcial porque varios valores permanecian dentro de las bandas, y esto debido a que las bandas resultaron muy amplias por la poca cantidad de registros. Finalmente, se logró proponer dos modelos, de los cuales los dos cumplieron los supuestos de los residuales, media cero, normalidad y blanquitud. El segundo modelo mostró menores valores en los criterios de información y así mismo obtuvo menores resultados en las medidas del Error Cuadratico Medio y el Error Porcentual Absoluto Medio por lo que se eligió este modelo para pronosticar el futuro, obteniendose buenos resultados.