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.