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.

2 comentarios:

  1. Este comentario ha sido eliminado por el autor.

    ResponderEliminar
  2. Hola Rosana.
    Mi nombre es Yonatan y soy de Perú.
    Te escribo porque estoy realizando mi tesis sobre Sistemas de Alerta Temprana con el modelo ARIMA usando imágenes MODIS.
    Me gustaría que me ayudes en cómo hacer un análisis de series de tiempo con ARIMA en R. Yo manejo R bastante bien, pero no sé como hacer mi analisis con mis datos.
    No sé si es posible por favor que me puedas brinda una ayuda ya sea por skype con la teoría de ARIMA y finalmente plantear el procedimiento en R con videos tutoriales u otro medio.

    Esta tesis es muy importante y te rogaría me contestaras pronto con una respuesta positiva o negativa.
    Gracias Rosana

    ResponderEliminar

Libros para descargar sobre Análisis de Series Temporales