#################################################
#################################################
#######  TEMA 3: MODELOS BOX-JENKINS  (1) #######
#################################################
#################################################


#######################################################
### Estudio pormenorizado de la serie9 (serie9.dat) ###
#######################################################

setwd("F:/German/Docencia/POP/08-09/Apuntes")

serie9 <- scan("serie9.dat")

max.ret <- trunc(length(serie9)/4)


# 1: IDENTIFICACIÓN (gráficos de la serie y de sus fas y fap)

  par(mfrow=c(3,1))
  plot.ts(serie9, xlab="Tiempo", ylab="")
  acf(serie9, xlab="Retardo", main="", lag.max=max.ret, xlim=c(2,max.ret), ylim=c(-1,1))
  pacf(serie9, xlab="Retardo", main="", lag.max=max.ret, xlim=c(2,max.ret), ylim=c(-1,1))

  # Observamos tendencia y componente estacional. Diferenciamos la serie regularmente:

  serie9.dif <- diff(serie9, lag=1)
  par(mfrow=c(3,1))
  plot.ts(serie9.dif, xlab="Tiempo", ylab="")
  acf(serie9.dif, xlab="Retardo", main="", lag.max=max.ret, xlim=c(2,max.ret), ylim=c(-1,1))
  pacf(serie9.dif, xlab="Retardo", main="", lag.max=max.ret, xlim=c(2,max.ret), ylim=c(-1,1))

  # Hemos eliminado la tendencia, aunque la componente estacional (s=12) persiste. Diferenciamos la serie estacionalmente:

  serie9.dif.dif12 <- diff(serie9.dif, lag=12)
  par(mfrow=c(3,1))
  plot.ts(serie9.dif.dif12, xlab="Tiempo", ylab="")
  acf(serie9.dif.dif12, xlab="Retardo", main="", lag.max=max.ret, xlim=c(2,max.ret), ylim=c(-1,1))
  pacf(serie9.dif.dif12, xlab="Retardo", main="", lag.max=max.ret, xlim=c(2,max.ret), ylim=c(-1,1))

  # La serie diferenciada (regular y estacionalmente) parece generada por un proceso estacionario

  # Sugerencia: ARIMA(0,1,1)x(0,1,1)_12 o ARIMA(1,1,0)x(0,1,1)_12 


# 2: ESTIMACIÓN

  # Si se ajusta/estima un modelo sin diferenciar (d=D=0), entonces por defecto la función arima() de R ajusta/estima un modelo con constante.
  # NOTA: En las salidas del programa, lo que identifica como "intercept" (esto es, como la estimación de la constante c del proceso)
  # es en realidad la estimación de la media del proceso (que coincide con la constante en el caso de los MA, pero no en el caso de 
  # procesos con componente AR)
  # Si deseamos ajustar/estimar un modelo sin diferenciar pero sin constante, debemos introducir en la funcion arima() el argumento include.mean=FALSE.

  # Si se ajusta/estima un modelo diferenciado (d o D no nulos) la función "arima" de R no permite introducir constante.
  # Puesto que la posibilidad de introducir constante no debe ser obviada, debemos realizar el ajuste/la estimación de un modo indirecto:
  # ajustamos/estimamos un ARMA con cte a la serie diferenciada. Si, como resultado, obtenemos que la media no es significativamente
  # distinta de cero, pasamos a ajustar un ARIMA (sin constante) a la serie sin diferenciar.

  # ARIMA(0,1,1)x(0,1,1)_12 

    # Ajustamos/estimamos un ARMA(0,1)x(0,1)_12 con cte a la serie diferenciada

    ma1.ma1.dif <- arima(serie9.dif.dif12, order=c(0,0,1), seasonal=list(order=c(0,0,1), period=12), optim.control=list(maxit=500)) 
    ma1.ma1.dif

    # Observamos que la media no es significativamente distinta de cero. 
    # Ajustamos entonces un ARIMA(0,1,1)x(0,1,1)_12 sin constante a la serie sin diferenciar

    arima011.011 <- arima(serie9, order=c(0,1,1), seasonal=list(order=c(0,1,1), period=12), optim.control=list(maxit=500)) 
    arima011.011


  # ARIMA(1,1,0)x(0,1,1)_12 

    # Ajustamos/estimamos un ARMA(1,0)x(0,1)_12 con cte a la serie diferenciada

    ar1.ma1.dif <- arima(serie9.dif.dif12, order=c(1,0,0), seasonal=list(order=c(0,0,1), period=12), optim.control=list(maxit=500)) 
    ar1.ma1.dif

    # Observamos que la media no es significativamente distinta de cero. 
    # Ajustamos entonces un ARIMA(1,1,0)x(0,1,1)_12 sin constante a la serie sin diferenciar

    arima110.011 <- arima(serie9, order=c(1,1,0), seasonal=list(order=c(0,1,1), period=12), optim.control=list(maxit=500)) 
    arima110.011


# 3: DIAGNOSIS

  # ARIMA(0,1,1)x(0,1,1)_12 

    # Gráficos

    par(mfrow=c(2,1))
    plot(arima011.011$residuals)
    qqnorm(arima011.011$residuals)

    # Contrastes

    tsdiag(arima011.011)
    t.test(arima011.011$residuals, mu=0)

    library(tseries)
    jarque.bera.test(arima011.011$residuals)
    shapiro.test(arima011.011$residuals)

  # Conclusión: El modelo  ARIMA(0,1,1)x(0,1,1)_12 sin constante resulta apropiado como generador de la serie, 
  # siendo además sus innovaciones gaussianas.


  # ARIMA(1,1,0)x(0,1,1)_12 

    # Gráficos

    par(mfrow=c(2,1))
    plot(arima110.011$residuals)
    qqnorm(arima110.011$residuals)

    # Contrastes
    tsdiag(arima110.011)
    t.test(arima110.011$residuals, mu=0)

    jarque.bera.test(arima110.011$residuals)
    shapiro.test(arima110.011$residuals)

  # Conclusión: El modelo  ARIMA(1,1,0)x(0,1,1)_12 sin constante resulta apropiado como generador de la serie, 
  # siendo además sus innovaciones Gaussianas.

  # COMPARACIÓN

  size <- length(serie9.dif.dif12)

    # (a) Criterio AIC

    AIC(arima011.011, arima110.011)

    # (b) Criterio AICc

    AIC(arima011.011, arima110.011, k=2*size/(size- 2 -2))

    # (c) Criterio BIC

    AIC(arima011.011, arima110.011, k=log(size))

  # Para cualquiera de los 3 criterios el menor valor es alcanzado por el ARIMA(1,1,0)x(0,1,1)_12. 
  # En cualquier caso, la diferencia entre sus valores es despreciable. 


# 4: IDENTIFICACIÓN (AIC y BIC)

  A <- array(0,c(3,3,3,3))
  B <- A

  for (p in 0:2) { for (q in 0:2) { for (P in 0:2) { for (Q in 0:2) { cat(p, q, P, Q, "\n"); a <- arima(serie9, order=c(p,1,q), seasonal=list(order=c(P,1,Q),   period=12), optim.control=list(maxit=500)); A[p+1,q+1,P+1,Q+1] <- a$aic; B[p+1,q+1,P+1,Q+1] <- AIC(a, k=log(size)) } } } }

  A <- A-min(A)
  B <- B-min(B)

  A

  # Modelo con menor AIC: ARIMA(1,1,1)x(0,1,1)_12 sin constante.
  # El modelo ARIMA(1,1,0)x(0,1,1)_12 sin constante tiene un AIC muy próximo al anterior (sus AICs distan menos de 2 unidades), y un parámetro menos. 
  # Este modelo ya lo habíamos chequeado.
  
  B

  # Modelo con menor BIC: ARIMA(1,1,0)x(0,1,1)_12 sin constante, el cual ya habíamos chequeado.
  # Además, no hay modelos con menos parámetros que el ARIMA(1,1,0)x(0,1,1)_12 cuyo BIC sea próximo al de éste (disten menos de 2 unidades).


# 5: PREDICCIONES BASADAS EN EL ARIMA(1,1,0)x(0,1,1)_12  AJUSTADO

  # Horizonte: 12 

    horizonte <- 12

    # Obtenemos las predicciones y sus desviaciones típicas

    serie9.pr <- predict(arima110.011, n.ahead=horizonte)
    serie9.pr

    # Construimos los intervalos de predicción (son válidos como consecuencia de la gaussianidad del ruido blanco del modelo)

    ext.inf <- serie9.pr$pred - 1.96*serie9.pr$se
    ext.sup <- serie9.pr$pred + 1.96*serie9.pr$se

    # Representamos graficamente los últimos 50 valores de la serie, junto con las predicciones y los correspondientes intervalos de predicción
    # hasta el horizonte establecido
  
    time.index <- (length(serie9)-49):length(serie9)
    min.t <- time.index[1]
    max.t <- max(time.index)+horizonte

    min.x <- min(serie9[time.index], ext.inf)
    max.x <- max(serie9[time.index], ext.sup)

    plot(time.index, serie9[time.index], type="o", xlim=c(min.t,max.t), ylim=c(min.x, max.x), xlab="Tiempo", ylab="")
    lines(serie9.pr$pred, col="blue", type="o")
    lines(ext.inf, col="red", lty="dashed")
    lines(ext.sup, col="red", lty="dashed")


  # Horizonte: 48

    horizonte <- 48

    serie9.pr <- predict(arima110.011, n.ahead=horizonte)
    serie9.pr

    ext.inf <- serie9.pr$pred - 1.96*serie9.pr$se
    ext.sup <- serie9.pr$pred + 1.96*serie9.pr$se

    time.index <- (length(serie9)-49):length(serie9)
    min.t <- time.index[1]
    max.t <- max(time.index)+horizonte

    min.x <- min(serie9[time.index], ext.inf)
    max.x <- max(serie9[time.index], ext.sup)

    plot(time.index, serie9[time.index], type="o", xlim=c(min.t,max.t), ylim=c(min.x, max.x), xlab="Tiempo", ylab="")
    lines(serie9.pr$pred, col="blue", type="o")
    lines(ext.inf, col="red", lty="dashed")
    lines(ext.sup, col="red", lty="dashed")