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


#####################################################
### Estudio pormenorizado de la serie turismo.txt ###
#####################################################

# turismo.txt: contiene la cantidad mensual de turistas que han visitado España el período enero 1995-agosto 2004

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

turismo <- scan("turismo.txt")


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

  par(mfrow=c(3,1))
  plot.ts(turismo, type="o", xlab="Mes", ylab="")
  acf(turismo, xlab="Retardo", main="", lag.max=length(turismo)/4, ylim=c(-1,1))
  pacf(turismo, xlab="Retardo", main="", lag.max=length(turismo)/4, ylim=c(-1,1))

  turismo.dif <- diff(turismo,1)
  par(mfrow=c(3,1))
  plot.ts(turismo.dif, type="o", xlab="Mes", ylab="")
  acf(turismo.dif, xlab="Retardo", main="", lag.max=length(turismo)/4, ylim=c(-1,1))
  pacf(turismo.dif, xlab="Retardo", main="", lag.max=length(turismo)/4, ylim=c(-1,1))

  turismo.dif.dif12 <- diff(turismo.dif,12)
  par(mfrow=c(3,1))
  plot.ts(turismo.dif.dif12, type="o", xlab="Mes", ylab="")
  acf(turismo.dif.dif12, xlab="Retardo", main="", lag.max=length(turismo)/4, ylim=c(-1,1))
  pacf(turismo.dif.dif12, xlab="Retardo", main="", lag.max=length(turismo)/4, ylim=c(-1,1))

  # SUGERENCIA: ARIMA(p,1,q)X(0,1,0)_12
  # NOTA: Quizás convenga darle estructura a la parte estacional, a la vista del valor fas(24), acompañado por un valor alto en la fas(25)
  # Ante tantas dudas, estudiamos modelos ARIMA(p,1,q)X(P,1,Q)_12


# 2: IDENTIFICACIÓN (BIC)

  size <- length(turismo.dif.dif12)

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



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

# ERROR 1 1 2 1


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

# ERROR 3 0 2 0


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


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

# ERROR 1 2 1 0


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

# ERROR 1 3 1 1


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


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

# ERROR 1 1 2 1


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


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


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


p <- 3
q <- 0
P <- 2
Q <- 1
a <- arima(turismo.dif.dif12, order=c(p,0,q), seasonal=list(order=c(P,0,Q), period=12), optim.control=list(maxit=500))
# ERROR

p <- 3
q <- 0
P <- 2
Q <- 2
a <- arima(turismo.dif.dif12, order=c(p,0,q), seasonal=list(order=c(P,0,Q), period=12), optim.control=list(maxit=500))
# ERROR



A[A==0] <- NaN

A <- A-min(A, na.rm=TRUE)
A


# Menor BIC: ARIMA(0,1,1)x(0,1,2)_12, seguido del ARIMA(2,1,3)x(0,1,1)_12 (a 0.91 unidades) y del ARIMA(0,1,1)x(2,1,2)_12 (a 1.58 unidades)
# Todos ellos con cte.


# El ARIMA(0,1,1)x(0,1,2)_12 con cte, aparte de tener el menor BIC, también es el más sencillo (menor cantidad de parámetros).
# Por tanto, si supera el análisis de residuos, sería el preferido entre los anteriores.

# 3: DIAGNOSIS

  # ARIMA(0,1,1)x(0,1,2)_12 con cte

    # Ajustamos el modelo

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


    # Gráficos

    par(mfrow=c(2,1))
    plot(ma1.ma2.dif12.dif.turismo$residuals)
    qqnorm(ma1.ma2.dif12.dif.turismo$residuals)

    # Contrastes
    tsdiag(ma1.ma2.dif12.dif.turismo, gof.lag=12)
    t.test(ma1.ma2.dif12.dif.turismo$residuals, mu=0)

    library(tseries)
    jarque.bera.test(ma1.ma2.dif12.dif.turismo$residuals)
    shapiro.test(ma1.ma2.dif12.dif.turismo$residuals)

    # Conclusión: El modelo  ARIMA(0,1,1)x(0,1,2)_12 con cte resulta apropiado como generador de la serie turismo, 
    # siendo sus innovaciones no gaussianas.


  
 # Puesto que uno de los órdenes del ARIMA(0,1,1)x(0,1,2)_12 coincide con el máximo valor que se le permitió (Q=2), 
   vamos a obtener el BIC de los modelos con Q=3.

 B <- array(0,c(4,4,3,1))


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

# ERROR 1 1 1 3


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

# ERROR 3 0 2 3


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

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


p <- 1
q <- 1
P <- 2
Q <- 3
a <- arima(turismo.dif.dif12, order=c(p,0,q), seasonal=list(order=c(P,0,Q), period=12), optim.control=list(maxit=500))
# ERROR


B[B==0] <- NaN

B

# Calculamos el BIC del mejor de los modelos ajustados en la primera etapa;
# esto es, el BIC del ARIMA(0,1,1)x(0,1,2)_12 con cte.

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

 AIC(a, k=log(size))

# Ninguno de los modelos propuestos con Q=3 supera al ARIMA(0,1,1)x(0,1,2)_12 con cte.

# Por tanto, éste es finalmente el seleccionado


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


  # Horizonte: 12 #

  horizonte <- 12


# Predicciones del proceso diferenciado

  dif12.dif.turismo.predicciones <- predict(ma1.ma2.dif12.dif.turismo, n.ahead=horizonte)$pred

  # Unimos la serie de observaciones diferenciadas con las predicciones obtenidas

  dif12.dif.turismo.y.predicciones <- c(turismo.dif.dif12, dif12.dif.turismo.predicciones)


# Deshacemos la diferencia estacional

  dif.turismo.y.predicciones <- diffinv( dif12.dif.turismo.y.predicciones, lag=12, xi=turismo.dif[1:12])


# Deshacemos la diferencia regular

  turismo.y.predicciones <- diffinv( dif.turismo.y.predicciones, lag=1, xi=turismo[1])



# Extraemos las predicciones

  pred.turismo<- turismo.y.predicciones[-(1:length(turismo))]

  

# Realizamos los gráficos

  min.t <- 1
  max.t <- length(turismo)+horizonte

  min.x <- min(turismo, pred.turismo)
  max.x <- max(turismo, pred.turismo)

  min.t.pred <- length(turismo)+1

  # Puesto que no tenemos normalidad, no disponemos de base teórica para la construcción de los intervalos de predicción.
  # Por tanto, no los construimos

  plot(turismo, type="o", xlim=c(min.t,max.t), ylim=c(min.x, max.x), xlab="Año", ylab="")
  lines(min.t.pred:max.t, pred.turismo, col="blue", type="o")