viernes, 7 de enero de 2011

Construcción de modelos ARIMA en R

#Construcción de modelos ARIMA
#1) graficamos los datos
gnp96 = read.table("c:/data/gnp96.dat") #Analizamos el producto nacional bruto (GNP) desde 1947(1)-2003(3), n=223
gnp = ts(gnp96[,2], start=1947, frequency=4) #construimos la serie
plot(gnp) #graficamos la serie
acf(gnp, 50) #graficamos el ACF y observamos que los datos no son estacionarios

#2) aplicamos posibles transformaciones a los datos
plot(diff(gnp)) #la diferenciación remueve la tendencia pero se observa aún algun tipo de tendencia remanente. La variabilidad en la segunda mitad es mayor que en la primer mitad de la serie.
gnpgr = diff(log(gnp)) # calculamos la tasa de crecimiento X(t)=diff(log(Y(t)))
plot.ts(gnpgr) # #parece ser un proceso estable

#3) identificamos los órdenes de dependencia del modelo
par(mfrow=c(2,1));acf(gnpgr, 24);pacf(gnpgr, 24) #Existen 2 interpretaciones posibles: 1) el ACF se corta en lag=2 y el PACf decrece, proponemos un modelo MA(2) para X(t) o ARIMA(0,1,2) para log(Y(t)) o 2) el ACF decrece y el PACF se corta en lag=1, identificándose un modelo AR(1) para X(t) o ARIMA(1,1,0) para log(Y(t)).

#4) estimamos los parámetros
# Como análisis preliminar ajustamos los dos modelos MA(2) y AR(1) por ML:
(gnpgr.ar = arima(gnpgr, order = c(1, 0, 0))) # potential problem here (see R Issue 1)
(gnpgr.ma = arima(gnpgr, order = c(0, 0, 2)))
ARMAtoMA(ar=.35, ma=0, 10) # devuelve los estimadores psi

#5) diagnosis: análisis de residuos y comparación de modelos
tsdiag(gnpgr.ar, gof.lag=20) #grafican los residuos estándar, el ACF de los residuos, el histograma y Q-Q plot de los residuales.
tsdiag(gnpgr.ma, gof.lag=20) #observamos que los residuos estándar no muestran un patrón obvio pero sí outliers, mientras que el ACF no muestra desvíos aparentes de los supuestos del modelo, y el Q-estadístico no es significativo. El histograma y Q-Q-plot de los residuos indica que estos son cercanos a la normalidad excepto algunos valores extremos en las colas.
hist(gnpgr.ma$resid, br=12); qqnorm(gnpgr.ma$resid); shapiro.test(gnpgr.ma$resid) #Sin embargo, el test de Shapiro-Wilks produce un p-valor de .003, indicando que los residuos no son normales
#el modelo MA(2) parece ajustarse bien a los datos salvo porque la distribución tiene colas que se apartan de la distribución normal

#6) elección del modelo, las técnicas más apropiada son el AIC, AICc y SIC.
n = length(gnpgr) # tamaño muestral
kar = length(gnpgr.ar$coef) # número de parámetros en el modelo AR(1)
sar=gnpgr.ar$sigma2 # estimación ML del sigma^2
kma = length(gnpgr.ma$coef) # número de parámetros en el modelo MA(2)
sma=gnpgr.ma$sigma2 # estimación ML del sigma^2

gnpgr.ar$aic #AIC, también se puede calcular como: log(sar) + (n+2*kar)/n # AR(1) -8.294
gnpgr.ma$aic #AIC. log(sma) + (n+2*kma)/n # MA(2) -8.298

log(sar) + (n+kar)/(n-kar-2) #AICc del AR(1) -8.285
log(sma) + (n+kma)/(n-kma-2) #AICc del MA(2) -8.288

log(sar) + kar*log(n)/n #BIO del AR(1) -9.264
log(sma) + kma*log(n)/n #BIC del MA(2) -9.252
#vemos que el AIC y AICc prefieren el modelo MA(2) mientras que el BIC selecciona el AR(1). Retenemos el modelo AR(1) porque es más simple (menor orden) y más fácil de trabajar.