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


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

# temperatura.txt: contiene la temperatura media mensual en la ciudad de La Coruña durante el período enero 1998-diciembre 2002

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

temperatura <- scan("temperatura.txt")

# Las 12 últimas observaciones no las utilizaremos en el ajuste. Posteriormente, serán utilizadas
# para comprobar la calidad de las predicciones.

temperatura1 <- temperatura[-((length(temperatura)-11):length(temperatura))]
temperatura2 <- temperatura[(length(temperatura)-11):length(temperatura)]


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

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

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

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


# 2: ESTIMACIÓN

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

  # La media del proceso diferenciado no resulta significativamente distinta de cero.  

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


# 3: DIAGNOSIS

  # Gráficos

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

  # Contrastes

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

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

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


# 4: IDENTIFICACIÓN (BIC)

  size <- length(temperatura1.dif12)

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

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

  A <- A-min(A)
  A


# Menor BIC: ARIMA(1,0,1)x(0,1,1)_12, seguido muy de cerca por el ARIMA(1,0,0)x(0,1,1)_12.

# Puesto que uno de los órdenes coincide con el máximo valor que se le permitió (q=1), vamos a obtener el BIC de los modelos con q=2.

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

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

B

# Calculamos el BIC del mejor de los modelos ajustados anteriormente: ARIMA(1,0,1)x(0,1,1)_12

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

AIC(a, k=log(size))


# Puesto que resulta mejor que cualquiera con q=2, y teniendo en cuenta la otra parte del estudio hecho,
# el modelo seleccionado es el ARIMA(1,0,0)x(0,1,1)_12 (es más sencillo que el del menor BIC: ARIMA(1,0,1)x(0,1,1)_12, y sus valores BIC son equivalentes).
# Este modelo ya había sido chequeado. 



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

  # Horizonte: 12 

  horizonte <- 12

  temperatura1.pr <- predict(arima100.011, n.ahead=horizonte)
  temperatura1.pr

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

  time.index <- 1:length(temperatura1)

  min.t <- time.index[1]
  max.t <- max(time.index)+horizonte

  min.x <- min(temperatura1, temperatura2, ext.inf)
  max.x <- max(temperatura1, temperatura2, ext.sup)

  plot(time.index, temperatura1, type="o", xlim=c(min.t,max.t), ylim=c(min.x, max.x), xlab="Mes", ylab="")
  lines(temperatura1.pr$pred, col="blue", type="o")
  lines(ext.inf, col="red", lty="dashed")
  lines(ext.sup, col="red", lty="dashed")
  lines((max.t-horizonte+1):max.t, temperatura2, col="green4", type="o")



####################################################
### Estudio pormenorizado de la serie tabaco.dat ###
####################################################

# tabaco.dat: contiene la producción anual de tabaco en EE.UU. durante el período  1872-1983

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

tabaco <- scan("tabaco.dat")


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

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

  plot.ts(tabaco, type="o", xlab="Año", ylab="")

  # El gráfico de la serie frente al tiempo sugiere que la variabilidad depende del nivel.
  # Por tanto, proponemos aplicarle a la serie una transformación de Box-Cox. 

  # El parámetro lambda de dicha transformación se elije de modo que la desviación típica sea una función potencial (potencia: 1-lambda) de la media.
  # Si bien a la vista del gráfico de secuencia parece que una transformación adecuada sería la logarítmica (lambda=0), podemos asesorarnos acerca
  # del valor de lambda a través de:

     # (a) El gráfico de la desviación típica estimada frente a la media estimada 
     #     (si observamos relación lineal tendríamos que 1-lambda=1 es apropiado)
     #     Estimamos la desviación típica y la media en cada tramo de, por ejemplo, 5 observaciones (si tuviésemos componente estacional de período s,      #      # utilizaríamos tramos de s observaciones). 
  
       aux <- matrix(tabaco,length(tabaco)/5, 5, byrow=TRUE)

       media <- apply(aux, 1, mean)
       des.tip <- apply(aux, 1, sd)

       plot(media, des.tip)

    # (b) Ajustamos una recta a los pares (ln(media), ln(desviación típica)) 
    # (un valor apropiado sería 1-lambda=pendiente de dicha recta)

      ln.media <- log(media)
      ln.des.tip <- log(des.tip)

      plot(ln.media, ln.des.tip)

      lm(ln.des.tip ~ ln.media)


  # Puesto que se aprecia cierta relación lineal (1-lambda=0.827 ~ 1), podemos considerar lambda=0 (esto es, la transformación logarítmica).

  ln.tabaco <- log(tabaco)

  plot.ts(ln.tabaco, type="o", xlab="Año", ylab="")

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

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

  # Sugerencia: ARIMA(0,1,1) o ARIMA(1,1,0) para ln(tabaco)


# 2: ESTIMACIÓN

  # ARIMA(0,1,1)

    ma1.dif.ln.tabaco <- arima(dif.ln.tabaco, order=c(0,0,1), optim.control=list(maxit=500)) 
    ma1.dif.ln.tabaco

    # La media del MA(1) resulta significativamente distinta de cero.  

  # ARIMA(1,1,0)

    ar1.dif.ln.tabaco <- arima(dif.ln.tabaco, order=c(1,0,0), optim.control=list(maxit=500)) 
    ar1.dif.ln.tabaco

    # La media del AR(1) no resulta significativamente distinta de cero.  

    arima110.ln.tabaco <- arima(ln.tabaco, order=c(1,1,0), optim.control=list(maxit=500)) 
    arima110.ln.tabaco


# 3: DIAGNOSIS

  # ARIMA(0,1,1)

    # Gráficos

    par(mfrow=c(2,1))
    plot(ma1.dif.ln.tabaco$residuals)
    qqnorm(ma1.dif.ln.tabaco$residuals)

    # Contrastes

    tsdiag(ma1.dif.ln.tabaco)
    t.test(ma1.dif.ln.tabaco$residuals, mu=0)

    library(tseries)
    jarque.bera.test(ma1.dif.ln.tabaco$residuals)
    shapiro.test(ma1.dif.ln.tabaco$residuals)

    # Conclusión: El modelo  ARIMA(0,1,1) con constante resulta apropiado como generador de la serie ln(tabaco), 
    # aunque sus innovaciones no son gaussianas.

  # ARIMA(1,1,0)

    # Gráficos

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

    # Contrastes
    tsdiag(arima110.ln.tabaco)
    t.test(arima110.ln.tabaco$residuals, mu=0)
    jarque.bera.test(arima110.ln.tabaco$residuals)
    shapiro.test(arima110.ln.tabaco$residuals)

    # Conclusión: El modelo  ARIMA(1,1,0) sin constante no supera claramente los contrastes de independencia. 
    # Por tanto, descartamos dicho modelo como generador de la serie ln(tabaco).


# 4: IDENTIFICACIÓN (BIC)

  size <- length(dif.ln.tabaco)

  A <- matrix(0,4,4)

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

  A <- A-min(A)
  A

  # Modelo con menor BIC: el ARIMA(0,1,1) con constante, el cual ya había sido chequeado. 
  # Además, los demás modelos son bastantes peores que éste.


# 5: PREDICCIONES BASADAS EN EL ARIMA(0,1,1) AJUSTADO (CON CONSTANTE)

  # Horizonte: 5 

  horizonte <- 5

  # Predicciones del proceso diferenciado

  dif.ln.tabaco.predicciones <- predict(ma1.dif.ln.tabaco, n.ahead=horizonte)$pred

  # Unimos la serie de observaciones (ln) diferenciadas con las predicciones obtenidas

  dif.ln.tabaco.y.predicciones <- c(dif.ln.tabaco, dif.ln.tabaco.predicciones)

  # Deshacemos la diferencia

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

  # Extraemos las predicciones

  pred.ln.tabaco <- ln.tabaco.y.pred[-(1:length(ln.tabaco))]

  # Deshacemos el ln

  pred.tabaco <- exp(pred.ln.tabaco)

  # Realizamos los gráficos

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

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

  min.t.pred <- length(tabaco)+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(tabaco, 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.tabaco, col="blue", type="o")