# El archivo "ibex35.ts.txt" contiene los rendimientos diarios (en porcentaje), excepto fines de semana y festivos, 
# del IBEX 35 desde el 13-11-2003 hasta el 11-11-2004

# Construiremos un modelo para su volatilidad. 
# Posteriormente, dicho modelo será utilizado para realizar predicciones de la volatilidad.

 setwd("F:/German/Docencia/POP/09-10/Alumnos/Volatilidad")

# setwd("H:/Datos_F_PC/German/Docencia/POP/09-10/Alumnos/Volatilidad")

# setwd("C:/MASTER-Profesores/Germán/09-10/Volatilidad")


##############################
### Serie "IBEX 35"
##############################

ibex35 <- scan("ibex35.txt")

ibex35.ts <- ts(ibex35)

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

  # Gráficos exploratorios
  
  windows()
  plot(ibex35.ts, xlab="Tiempo", ylab="IBEX 35")

  windows()
  qqnorm(ibex35.ts)
  qqline(ibex35.ts)

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

  # Fas y fap de la serie al cuadrado  
  windows()
  par(mfrow=c(2,1))
  acf(ibex35.ts^2, xlab="Retardo", main="", lag.max=trunc(length(ibex35.ts)/8), xlim=c(1,trunc(length(ibex35.ts)/8)), ylim=c(-1,1))
  pacf(ibex35.ts^2, xlab="Retardo", main="", lag.max=trunc(length(ibex35.ts)/8), xlim=c(1,trunc(length(ibex35.ts)/8)), ylim=c(-1,1))

  # Fas y fap del valor absoluto de la serie  
  windows()
  par(mfrow=c(2,1))
  acf(abs(ibex35.ts), xlab="Retardo", main="", lag.max=trunc(length(ibex35.ts)/8), xlim=c(1,trunc(length(ibex35.ts)/8)), ylim=c(-1,1))
  pacf(abs(ibex35.ts), xlab="Retardo", main="", lag.max=trunc(length(ibex35.ts)/8), xlim=c(1,trunc(length(ibex35.ts)/8)), ylim=c(-1,1))

# 2: BUSCAMOS EL MEJOR ARMA (BIC) COMO GENERADOR DE LA SERIE: 
#    el estudio anterior sugiere incorrelación pero no independencia, y colas pesadas.
#    El GARCH es ruido blanco: esto es, un ARMA(0,0) sin cte.

   source("F:/German/Docencia/POP/09-10/Alumnos/Box_Jenkins/Box_Jenkins_Seleccion.R")

   mejor.arma(x=ibex35.ts, ord.max=c(5,5,0,0), incluir.media="SI", metodo=NULL, criterio="BIC", dist.max.crit=2)

# Modelo tentativo (BIC, máxima verosimilitud): ARMA(0,0)

# 3: AJUSTE

  ajuste.arma <- arima(ibex35.ts, order=c(0,0,0), include.mean=TRUE, optim.control=list(maxit=500))
  ajuste.arma

# Cte no significativamente distinta de cero

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

# 4. ANÁLISIS DE RESIDUOS (nótese que los residuos coinciden con la serie; esto es, residuals(ajuste.arma)=ibex35.ts)

  tsdiag(ajuste.arma)
  t.test(residuals(ajuste.arma), mu=0)

  library(tseries)
  jarque.bera.test(residuals(ajuste.arma))
  shapiro.test(residuals(ajuste.arma))

  # Conclusión: Se puede considerar que la serie ibex35.ts 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.


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

############ Utilizamos el criterio BIC para proponer valores para r y s
############ La función "mejor.garch" (que se encuentra en el script "Volatilidad_Seleccion.R") lleva a cabo este proceso

  source("Volatilidad_Seleccion.R")

  mejor.garch(x=ibex35.ts, ord.max=c(4,4), criterio="BIC", dist.max.crit=2)

# Modelo tentativo: GARCH(3,0) o, lo que es lo mismo, ARCH(3)

 
# 6: ESTIMACIÓN DEL GARCH IDENTIFICADO

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

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

  # Desviaciones típicas de los estimadores

   sqrt(diag(arch3$vcov))

  # Coeficientes no significativos: el alpha1 y alpha3

  arch3$coef / (1.96*sqrt(diag(arch3$vcov)))

  # Los parámetros alpha_1 y alpha_3 resultan no significativamente distintos de cero.
  # La función garch de R no permite fijar a cero coeficientes diferentes del último. 
  # Lo que haremos es mantener el alpha1 y eliminar el alpha3:
  # Pasamos entonces a estimar un ARCH(2) (es el modelo que resulta de fijar alpha3=0)

  arch2 <- garch(ibex35.ts, order=c(0,2), trace=FALSE)

  arch2$coef

  sqrt(diag(arch2$vcov))
 
  arch2$coef / (1.96*sqrt(diag(arch2$vcov)))

  # Aunque el parámetro alpha_1 sigue siendo no significativamente distinto de cero, lo mantenemos (debido a las carencias de la función "garch").


# 7: 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

  windows()
  par(mfrow=c(2,1))
  plot(residuals(arch2), xlab="Tiempo", ylab="")
  abline(h=0)

  qqnorm(residuals(arch2))
  qqline(residuals(arch2))

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

       # (b) Media 0

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

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

       # (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) -2
           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.
  # Su ruido Z_t no es gaussiano



# 8: 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.ts, 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.ts)+1, dt.cond.y.pred[,1][length(dt.cond.y.pred[,1])], type="p", pch=19, col="red")













################################################
# ETAPA: IDENTIFICACIÓN de los órdenes (p, d, q)
#################################################

# Gráfico de secuencia
windows()
plot(nilo.ts, xlab="Año", ylab="Nivel del Nilo")

# Fas y fap estimadas
windows()
par(mfrow=c(2,1))
acf(nilo.ts, lag.max=75, main="")
pacf(nilo.ts, lag.max=75, main="")

windows()
plot(log(1:75), log(acf(nilo.ts, lag.max=75, plot=FALSE)$acf[-1]), ylab="log(fas)", xlab="log(retardo)")

# Modelos tentativo: FARIMA(p,d,q)

############ Utilizamos el criterio BIC para proponer valores para p, d y q
############ La función "mejor.farima" (que se encuentra en el script "Box_Jenkins_SeleccionFarima.R") lleva a cabo este proceso

source("Memoria_Larga_Seleccion.R")

mejor.farima(x=nilo.ts, ord.max=c(2,2), criterio="BIC", dist.max.crit=2)

# Modelo tentativo (BIC): FARIMA(0,0.3933,0)

#############################
# ETAPA: ESTIMACIÓN O AJUSTE 
#############################

### Ajuste

source("fracdiff.master.R")

ajuste <- fracdiff.master(x=nilo.ts)

ajuste$coef

ajuste$intconf

# El parámetro d es significativamente distinto de cero

mean(nilo.ts)

##############################
# ETAPA: ANÁLISIS DE RESIDUOS
##############################

residuos <- ajuste$residuals

# Gráfico de secuencia de los residuos del ajuste arma anterior
windows()
plot(residuos, type="o", xlab="Año", ylab="Residuos")
abline(h=0)

# Gráfico Q-Q normal
windows()
qqnorm(residuos)
qqline(residuos)

# Fas estimadas y Ljung-Box
windows()
tsdiag(ajuste$ajuste.arma, gof.lag=20)

# mu_a = 0
t.test(residuos, mu=0)

# Normalidad
library(tseries)
jarque.bera.test(residuos)
shapiro.test(residuos)

## Conclusión: El modelo ajustado FARIMA(0,0.3933,0)
 con media mean(nilo.ts)=1148.125 puede ser utilizado como generador de la serie "nilo". Las innovaciones NO son gaussianas.


####################
# ETAPA: PREDICCIÓN
####################

source("predict.farima.R")

nilo.pr <- predict.farima(x=nilo.ts, d=ajuste$d, horizonte=50)

windows()
ts.plot(cbind(window(nilo.ts, start=614), nilo.pr$pred), col=c("black", "blue"), type="o", xlab="Año", ylab="Nivel")
abline(h=mean(nilo.ts), col="green")
legend(x=614, y=max(window(nilo.ts, start=614)), legend = c("Serie temporal", "Predicciones", "Media del FARIMA"), col=c("black", "blue", "green"), lty=1, pch=c(1,1,NA), bty="n")



#########################################################################################################
#
# A continuación, modelizamos la misma serie (nilo) a través de un proceso de memoria corta (ARMA).
# Después, compararemos las preducciones obtenidas a través del FARIMA con las correspondientes al ARMA.
#
#########################################################################################################

############ Utilizamos el criterio BIC para construir modelos tentativos
############ La función "mejor.arma" (que se encuentra en el script "Box_Jenkins_Seleccion.R") lleva a cabo este proceso

source("F:/German/Docencia/POP/09-10/Alumnos/Box_Jenkins/Box_Jenkins_Seleccion.R")

mejor.arma(x=nilo.ts, ord.max=c(8,8,0,0), incluir.media="SI", metodo=NULL, criterio="BIC", dist.max.crit=2)

# Modelo tentativo (BIC, máxima verosimilitud): ARIMA(2,0,1)

#############################
# ETAPA: ESTIMACIÓN O AJUSTE 
#############################

### Ajuste por máxima verosimilitud: ARIMA(2,0,1)

source("F:/German/Docencia/POP/09-10/Alumnos/Box_Jenkins/Box_Jenkins_Coef_Significativos.R")

ajuste.arma21 <- arma.signif(x=nilo.ts, orden=c(2,1,0,0), metodo=NULL)
ajuste.arma21

# Modelo tentativo BIC: ARIMA(2,0,1) CON constante

##############################
# ETAPA: ANÁLISIS DE RESIDUOS
##############################

# Gráfico de secuencia de los residuos del ajuste arma anterior
windows()
plot(residuals(ajuste.arma21), type="o", xlab="Año", ylab="Residuos")
abline(h=0)

# Gráfico Q-Q normal
windows()
qqnorm(residuals(ajuste.arma21))
qqline(residuals(ajuste.arma21))

# Fas estimadas y Ljung-Box
windows()
tsdiag(ajuste.arma21, gof.lag=20)

# mu_a = 0
t.test(residuals(ajuste.arma21), mu=0)

# Normalidad
library(tseries)
jarque.bera.test(residuals(ajuste.arma21))
shapiro.test(residuals(ajuste.arma21))

## Conclusión: El modelo ajustado ARIMA(2,0,1) CON cte 
##             puede ser utilizado como generador de la serie "nilo". Las innovaciones NO son gaussianas.


####################
# ETAPA: PREDICCIÓN
####################

############ Utilizando el modelo ajustado ARIMA(2,0,1) CON cte, predecimos los próximos 50 valores del proceso

# Aunque el modelo contiene constante, no se ha diferenciado. Por tanto, podemos utilizar la función "predict"

# Predicción de futuros valores del nilo.ts

nilo.arma.pr <- predict(ajuste.arma21, n.ahead=50)

media <- ajuste.arma21$coef[length(ajuste.arma21$coef)]

windows()
ts.plot(cbind(window(nilo.ts, start=614), nilo.arma.pr$pred), col=c("black", "red"), type="o", xlab="Año", ylab="Nivel")
abline(h=media, col="green")
legend(x=614, y=max(window(nilo.ts, start=614)), legend = c("Serie temporal", "Predicciones", "Media del ARMA"), col=c("black", "red", "green"), lty=1, pch=c(1,1,NA), bty="n")


##################################################################################
# COMPARACIÓN GRÁFICA DE LAS PREDICCIONES OBTENIDAS POR LOS MODELOS FARIMA Y ARMA
##################################################################################

windows()
par(mfrow=c(1,2))

ts.plot(cbind(window(nilo.ts, start=614), nilo.pr$pred), col=c("black", "blue"), type="o", xlab="Año", ylab="Nivel")
abline(h=mean(nilo.ts), col="green")
legend(x=614, y=max(window(nilo.ts, start=614)), legend = c("Serie temporal", "Predicciones", "Media del FARIMA"), col=c("black", "blue", "green"), lty=1, pch=c(1,1,NA), bty="n")

ts.plot(cbind(window(nilo.ts, start=614), nilo.arma.pr$pred), col=c("black", "red"), type="o", xlab="Año", ylab="Nivel")
abline(h=media, col="green")
legend(x=614, y=max(window(nilo.ts, start=614)), legend = c("Serie temporal", "Predicciones", "Media del ARMA"), col=c("black", "red", "green"), lty=1, pch=c(1,1,NA), bty="n")