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


#########################################################
### 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")


setwd("C:/MASTER-Profesores/Germán/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


#####################################################
### 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")

setwd("C:/MASTER-Profesores/Germán/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


####################################################
### 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")


setwd("C:/MASTER-Profesores/Germán/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 la pendiente de la recta ajustada es próxima a 1 (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)


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

# consumo_electricidad.txt: contiene el consumo mensual de electricidad en EE.UU. durante el periodo enero 1972-diciembre 2004


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


setwd("C:/MASTER-Profesores/Germán/08-09/Apuntes")


consumo <- scan("consumo_electricidad.txt")

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

consumo1 <- consumo[1:(length(consumo)-12)]
consumo2 <- consumo[-(1:(length(consumo)-12))]

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

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

  plot.ts(consumo1, type="o", xlab="Mes", 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 12 observaciones (observamos componente estacional con s=12). 
  
       aux <- matrix(consumo1,length(consumo1)/12, 12, 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 la pendiente de la recta ajustada es próxima a 1 (1-lambda=1.298 ~ 1), podemos considerar lambda=0 (esto es, la transformación logarítmica).


  ln.consumo1 <- log(consumo1)

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


  dif.ln.consumo1 <- diff(ln.consumo1, lag=1)

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


  dif12.dif.ln.consumo1 <- diff(dif.ln.consumo1, lag=12)

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


  # Sugerencia: ARIMA(0,1,2)x(0,1,1)_12 para ln(consumo1)


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

# viviendas.dat: contiene la cantidad mensual de viviendas iniciadas en Galicia durante el período enero 1987-diciembre 2000


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


setwd("C:/MASTER-Profesores/Germán/08-09/Apuntes")



viviendas <- scan("viviendas.dat")


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

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

  plot.ts(viviendas, type="o", xlab="Mes", ylab="")

  # Aunque no está demasiado claro, quizás la variabilidad dependa del nivel
  # (la zona de mayor variabilidad se encuentra al final de la serie, que coincide con la zona de mayor nivel)

  # Trataremos entonces de construir una transformación razonable (objetiva).
  # Para ello, ajustamos una recta a los pares (ln(media), ln(desviación típica)) 
  # (un valor apropiado sería 1-lambda=pendiente de dicha recta)

  aux <- matrix(viviendas,length(viviendas)/5, 5, byrow=TRUE)

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

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

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

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

  # Se tiene que 1-lambda=0.5495. Podríamos entonces considerar lambda=0.5 (esto es, transformar la serie a través de la raíz cuadrada).

  # El coeficiente de determinación (R^2) correspondiente al ajuste lineal que se acaba de realizar ha resultado ser 0.1403. 
  # Esto indica falta de linealidad, y por tanto no es esperable que la transformación de Box-Cox con el parámetro lambda propuesto funcione bien
  # (la transformación de Box-Cox funciona bien si hay relación potencial entre la desviación típica y la media o, lo que es lo mismo, 
  # si hay relación lineal entre los logaritmos de dichas medidas).
  # Por otra parte, se puede comprobar que, en este ejemplo concreto, pequeños cambios en la cantidad de datos utilizados para realizar la estimación de
  # la media y la desviación típica redundan en grandes cambios en el valor de lambda (esto es, el lambda depende enormemente de la cantidad de
  # observaciones utilizadas en las estimaciones). 
  # Esto puede considerarse un mal augurio en lo que a la fiabilidad de la transformación de Box-Cox se  refiere (siempre para este ejemplo concreto).

  # En cualquier caso, continuamos con el estudio.

  raiz.viviendas <- sqrt(viviendas)

  # La transformación de Box-Cox propone (x^lambda - 1)/lambda para estabilizar la varianza.
  # Obviamente, si la anterior transformación consigue homocedasticidad, también la consigue la transformación x^lambda.
  # Las constantes que acompañan a x^lambda en la transformación de Box-Cox se introducen únicamente para que al tender 
  # lambda a cero la transformación congerja a la transformación logarítmica.

  windows()
  par(mfrow=c(2,1))
  plot.ts(viviendas, type="o", xlab="Mes", ylab="viviendas")
  plot.ts(raiz.viviendas, type="o", xlab="Mes", ylab="raiz(viviendas)")

  # Recordemos los gráficos de las funciones raíz cuadrada y logaritmo neperiano:

  x <- seq(0.1,20,length=200)

  y.min <- min(sqrt(x), log(x))
  y.max <- max(sqrt(x), log(x))

  windows()
  plot(x, sqrt(x), ylim=c(y.min, y.max), type="l", col="blue", ylab="")
  lines(x, log(x), type="l", col="red", ylab="")
  abline(0,1)
  legend(10, 0, lty=1, c("x", "sqrt(x)", "ln(x)"), col=c("black", "blue", "red"))

  # Observemos cuál es el efecto de aplicar cada una de estas dos transformaciones:
  # Ambos transformaciones separan valores muy próximos a cero, y juntan valores alejados del cero:
  # Cuanto más próximos están de cero los valores, más los separan; cuanto más alejados están de cero los valores, más los juntan.
 
  # La diferencia principal entre ambas transformaciones radica en que la función logarítmica da lugar a una separación mayor 
  # (de los valores muy próximos a cero) que la raíz cuadrada, 
  # y también a un acercamiento mayor (de los valores alejados del cero) que la raíz cuadrada. 

  # Veamos ahora cómo estas funciones transforman nuestra serie:

  ln.viviendas <- log(viviendas)

  windows()
  par(mfrow=c(3,1))
  plot.ts(viviendas, type="o", xlab="Mes", ylab="viviendas")
  plot.ts(raiz.viviendas, type="o", xlab="Mes", ylab="raiz(viviendas)")
  plot.ts(ln.viviendas, type="o", xlab="Mes", ylab="ln(viviendas)")

  # Parece que la transformación logarítmica resulta más apropiada que la raíz cuadrada (aunque la transformación raíz cuadrada era la propuesta por
  # el ajuste lineal hecho anteriormente, arriba hemos indicado su poca fiabilidad en este ejemplo concreto).
  # La superioridad de la transformación logarítmica puede justificarse recordando que la logarítmica junta más los valores grandes 
  # de lo que lo hace la raíz cuadrada; además, en nuestra serie todos los valores están alejados del cero, 
  # estando el cambio de variabilidad en los más altos.

  # Otro punto a favor de la transformación logarítmica está en que, si después hay que diferenciar la serie transformada, se tiene que:

  # ln(x_t) - ln(x_(t-1)) = ln(x_t / x_(t-1)) = ln (1 + R_t)

  # donde

  # R_t = (x_t - x_(t-1)) / x_(t-1).

  # Entonces, si R_t es próximo a cero, se tiene que

  # ln(x_t) - ln(x_(t-1)) = ln (1 + R_t) ~ R_t

  # En muchas series financieras R_t se conoce como el rendimiento, y representa el cambio que experimenta
  # la serie original desde el instante t-1 hasta el instante t, con respecto a su valor en el instante t-1.
  # Por tanto, trabajar en unidades logarítmicas resulta prácticamente equivalente a trabajar con la serie de rendimientos.

  # Resumiendo, aunque al realizar el ajuste lineal la transformación propuesta resultó ser la función raíz cuadrada,
  # tenemos motivos para quedarnos con la transformación logarítmica.

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

  # El gráfico de las correlaciones simples muestrales muestra valores positivos que decrecen lentamente a cero (decrecimiento lineal).
  # Si sus valores fuesen próximos a 1, diferenciaríamos regularmente la serie (presencia de tendencia).
  # Si bien en el caso que nos ocupa dichos valores no son próximos a 1, el hecho de que sus valores sean positivos y presenten 
  # decrecimiento lento de forma lineal es suficiente para indicar no estacionariedad (y necesidad de diferenciar regularmente). 

  dif.ln.viviendas <- diff(ln.viviendas,1)

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

  # Sugerencia: ARIMA(0,1,1)x(1,0,0)_12 para el ln de la cantidad mensual de viviendas...



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

# colgate.dat: contiene la cuota de mercado semanal del dentífrico "Colgate" durante el período 01/enero/1958-30/abril/1963


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


setwd("C:/MASTER-Profesores/Germán/08-09/Apuntes")



colgate <- scan("colgate.dat")


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

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

# Las correlaciones simples estimadas de esta serie, si bien son positivas, presentan un decrecimiento todavía más lento que el decrecimiento lineal.
# Podemos decir que no acaban de desaparecer, sino que persisten en el tiempo (aunque con valores bajos). 
# Debe notarse que esta característica es nueva con respecto a lo visto hasta ahora. Series presentando esta propiedad serán estudiadas en el próximo tema.