###########################################################
#######     TEMA 5: MODELOS PARA LA VOLATILIDAD    ########
###########################################################

###############################################################
## Estudio pormenorizado de la serie del Ibex35 (ibex35.txt) ##
###############################################################


# ibex.txt: contiene los rendimientos diarios del Ibex 35 desde el 13-11-2003 hasta el 11-11-2004 (excepto fines de semana y festivos)

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


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


ibex35 <- scan("ibex35.txt")


# 1: DETECCIÓN DE LA PRESENCIA DE HETEROCEDASTICIDAD CONDICIONAL

  # Gráficos secuencial y Q-Q, en los que se observan características presentes en series generadas por procesos GARCH
  
  par(mfrow=c(1,2))
  plot.ts(ibex35, xlab="Tiempo", ylab="IBEX 35")
  qqnorm(ibex35)


  # Fas y fap muestrales de la serie y sus cuadrados.
  # En ellas se observan características presentes en series generadas por procesos GARCH
  
  par(mfrow=c(2,1))
  acf(ibex35, xlab="Retardo", main="", lag.max=trunc(length(ibex35)/8), xlim=c(1,trunc(length(ibex35)/8)), ylim=c(-1,1))
  pacf(ibex35, xlab="Retardo", main="", lag.max=trunc(length(ibex35)/8), xlim=c(1,trunc(length(ibex35)/8)), ylim=c(-1,1))

  par(mfrow=c(2,1))
  acf(ibex35^2, xlab="Retardo", main="", lag.max=trunc(length(ibex35)/8), xlim=c(1,trunc(length(ibex35)/8)), ylim=c(-1,1))
  pacf(ibex35^2, xlab="Retardo", main="", lag.max=trunc(length(ibex35)/8), xlim=c(1,trunc(length(ibex35)/8)), ylim=c(-1,1))

 
  # LA SERIE EN ESTUDIO... ¿PROVIENE DE UN PROCESO DE RUIDO BLANCO (EL PROSESO GARCH LO ES)?

    # La función tsdiag sólo se puede aplicar sobre objetos generados por un ajuste.
    # Por tanto, ajustamos un ARMA(0,0) sin constante a ibex35, y le aplicamos tsdiag a dicho ajuste
    # Nótese que en dicho ajuste únicamente se estima la varianza, y los residuos coinciden con la serie original ibex35.

  ruido.blanco.ibex35 <- arima(ibex35, order=c(0,0,0), include.mean=FALSE, optim.control=list(maxit=500))
  ruido.blanco.ibex35

  tsdiag(ruido.blanco.ibex35)
  t.test(ibex35, mu=0)

  library(tseries)
  jarque.bera.test(ibex35)
  shapiro.test(ibex35)

  # Conclusión: Se puede considerar que la serie ibex35 ha sido generada por un proceso de ruido blanco no gaussiano.
  # Recuérdese que bajo normalidad, incorrelación e independencia son conceptos equivalentes. Por tanto, puesto que el
  # GARCH es ruido blanco con estructura de dependencia, no puede ser un proceso Gaussiano.


  # ¿EXISTE ALGÚN ARMA DISTINTO DEL RUIDO BLANCO QUE RESULTE MÁS APROPIADO QUE ÉSTE PARA HABERLA GENERADO?

  size <- length(ibex35)

  A <- matrix(0,4,4)

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

  # Respuesta: NO


# 2: IDENTIFICACIÓN (BIC) DE UN PROCESO GARCH COMO GENERADOR DE LA SERIE ibex35.

  # La función "garch" (incluida en el paquete "tseries") de R utiliza los órdenes del GARCH en orden inverso al usual.

  size <- length(ibex35)
  A <- matrix(0,5,5)
  for (r in 0:4) { for (s in 0:4) { if (!((r==0) & (s==0)))  { A[r+1,s+1] <- AIC(garch(ibex35, order=c(s,r), trace=FALSE), k=log(size-max(r,s)))}}}

  A[1,1] <- NaN

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

  # El ARCH(3) es el que alcanza un menor BIC. Pasamos a estimarlo


# 3. ESTIMACIÓN DEL GARCH IDENTIFICADO

  arch3 <- garch(ibex35, order=c(0,3), trace=FALSE)
  arch3$coef

  # Desviaciones típicas de los estimadores

  sqrt(arch3$vcov[1,1])
  sqrt(arch3$vcov[2,2])
  sqrt(arch3$vcov[3,3])
  sqrt(arch3$vcov[4,4])

  # Los parámetros alpha_1 y alpha_3 resultan no significativamente distintos de cero.
  # Pasamos entonces a estimar un ARCH(2) (es el modelo que resulta de fijar alpha3=0)


  arch2 <- garch(ibex35, order=c(0,2), trace=FALSE)
  arch2$coef

  sqrt(arch2$vcov[1,1])
  sqrt(arch2$vcov[2,2])
  sqrt(arch2$vcov[3,3])

  # Aunque el parámetro alpha_1 sigue siendo no significativamente distinto de cero, lo mantenemos, pues la función "garch" no permite fijarlo a cero.


# 4: ANÁLISIS O CHEQUEO DE LOS RESIDUOS DEL ARCH(2)

  # A continuación, mostramos un análisis de residuos generado automáticamente por las 2 órdenes siguientes:

  summary(arch2)

  plot(arch2)


  # Pasamos ahora a completar el análisis de los residuos (y de sus cuadrados):

    # RESIDUOS
 
      # (a) Independencia

      # La función tsdiag sólo se puede aplicar sobre objetos generados por un ajuste (ajuste arima,..., pero no ajuste garch).
      # Por tanto, ajustamos un ARMA(0,0) sin constante a los residuos que queremos analizar, y le aplicamos tsdiag a dicho ajuste.
      
      aux <- arima(arch2$residuals[-(1:2)], order=c(0,0,0), include.mean=FALSE, optim.control=list(maxit=500))
      tsdiag(aux)

      # Otra forma de conseguir lo mismo, sería a através del siguiente código:

          # Obtenemos los p-valores del contraste de Ljung-Box
            H <- 10
            p.valores <- 0  
    
            for (i in 1:H) p.valores[i] <- Box.test(arch2$residuals[-(1:2)], lag=i, type="Ljung-Box")$p.value
            p.valores

          # Hacemos los gráficos
 
            par(mfrow=c(3,1))
            plot(arch2$residuals, type="l", xlab="Time", ylab="", main="")
            acf(arch2$residuals, main="", lag.max=H, na.action=na.pass)
            plot(p.valores, xlab="Lag", ylim=c(0,1))
            abline(0.05, 0, lty=2)
    

       # (b) Media 0

         t.test(arch2$residuals, mu=0)

       # (c) Normalidad
  
        jarque.bera.test(arch2$residuals[-(1:3)])
        shapiro.test(arch2$residuals)

       # (d) Varianza 1 (sólo es válido si los residuos son normales)
  
         # Construimos un IC al 95% para la varianza

           size.residuals <- length(arch2$residuals) -3
           numerador <- (size.residuals -1)*var(arch2$residuals, na.rm=TRUE)
           ext.inf <- numerador/qchisq(0.975, size.residuals -1)
           ext.sup <- numerador/qchisq(0.025, size.residuals -1)
           cat("(", ext.inf, " , ", ext.sup, ")", "\n")

         # Obtenemos el p-valor
   
           estadistico <- numerador
           p.value <- 2*min(pchisq(numerador,size.residuals -1), 1-pchisq(numerador,size.residuals -1))
           p.value


   # RESIDUOS^2

      # Independencia

      aux <- arima(arch2$residuals[-(1:2)]^2, order=c(0,0,0), include.mean=FALSE, optim.control=list(maxit=500))
      tsdiag(aux)

 
  # Conclusión: Un modelo ARCH(2) puede ser considerado apropiado como generador de la serie de rendimientos del IBEX 35.
  # Deberíamos ser un tanto escépticos ante la hipótesis de ruido gaussiano


# 5: SERIE IBEX35 vs RESIDUOS ARCH(2) 

  # A continuación, mostramos algunos gráficos comparando el comportamiento de la serie de rendimientos con la de los residuos.

    # Gráficos secuenciales

     y.min <- min(ibex35, arch2$residuals, na.rm=TRUE)
     y.max <- max(ibex35, arch2$residuals, na.rm=TRUE)

     plot.ts(ibex35, ylab="", ylim=c(y.min,y.max))
     plot.ts(arch2$residuals, ylab="", ylim=c(y.min,y.max))


  # Q-Q normal

     qqnorm(ibex35, ylim=c(y.min,y.max), main="")
     qqnorm(arch2$residuals, ylim=c(y.min,y.max), main="")


  # fas de los cuadrados

    acf(ibex35^2, xlab="Retardo", main="", lag.max=30, xlim=c(1,30), ylim=c(-0.2,1))
    acf(arch2$residuals^2, xlab="Retardo", lag.max=30, xlim=c(1,30), ylim=c(-0.2,1), na.action=na.pass)

# 6: PREDICCIÓN

  # Comparamos en primer lugar la serie de rendimientos del IBEX 35 con su desviación típica condicional ajustada

   par(mfrow=c(2,1))
   plot.ts(ibex35, ylab="IBEX 35")
   plot.ts(arch2$fitted.values[,1], ylab="Des. típ. cond")

  # Por último, predecimos a horizonte 1 la desviación típica condicional

   dt.cond.y.pred <- predict(arch2, genuine=TRUE) 
   plot.ts(dt.cond.y.pred[,1][-length(dt.cond.y.pred[,1])], ylab="Des. típ. cond")
   points(length(ibex35)+1, dt.cond.y.pred[,1][length(dt.cond.y.pred[,1])], type="p", pch=19, col="red")