#R : R 3.4.1
# Ex1_gspc_garch_R.txt

# source("Z:/LV/QF_FinancialEconometrics/CHP1_GARCH_DF_Endo_Caus_IV/EXERCISES/Ex1_gspc_garch_R.txt") 
#setwd("C:/MH/LV/QF_FinancialEconometrics/CHP1_GARCH_DF_Endo_Caus_IV/Daten/")
#setwd("Z:/LV/QF_FinancialEconometrics/CHP1_GARCH_DF_Endo_Caus_IV/Daten/")

library("tseries")
library("zoo")
library("rugarch")
library("urca")

source("BasicStatistics_R.txt")

# FX: USD/EUR
tab <- read.table(file="fx_usd_eur.txt", header = FALSE, sep = "", 
           dec = ".", col.names = c("day", "fx" ) )
n <- length(tab[,1]) # 754, ordering from newest to oldest
day <- rev(tab[,1])
fx  <- rev(tab[,2])
fx.Date <- as.Date("2014-09-01") + c(1:n) - 1

#write.table(fx, file="nix.txt", row.names=FALSE, col.names=FALSE )

fx.z        <- as.zoo(fx)
index(fx.z) <- fx.Date
d_fx.z      <- diff(fx.z)
d_fx2.z     <- ( d_fx.z - mean(d_fx.z) )^2

y   <- fx.z
dy  <- d_fx.z
dy2 <- d_fx2.z

### Plotting univariate data,      Press ESC to continue
plot(y,   main = "y");          Sys.sleep(5)   # exchange rate, oanda 
plot(dy,  main = "Delta y");    Sys.sleep(5)
plot(dy2, main = "(Delat y)^2" )  


### Testing for white noise 
acf(y);                     # acf: very slow decay nonstat?
acf(dy);                    # acf: white noise ?
 
LjungBox(y,10)   # Ljung-Box tests
LjungBox(dy,10) 

##Dickey-Fuller Unit root test: level and 1st difference
# tau* test for  H0:                (1-rho) = 0   in various models
# phi1  test for H0: const=0   and  (1-rho) = 0   in model w const/drift
# phi3  test for H0: no trend  and  (1-rho) = 0   in trend model
## type = "none" or "drift" or "trend"
## selectlags= "Fixed" or "AIC", ... 
## relevant statistic: tau

# y
summary( ur.df(y,  type="none",  lags=5, selectlags="AIC")  )
summary( ur.df(y,  type="drift", lags=5, selectlags="AIC")  )
summary( ur.df(y,  type="trend", lags=5, selectlags="AIC")  )

# dy
summary( ur.df(dy,  type="none",  lags=5, selectlags="AIC")  )
summary( ur.df(dy,  type="drift", lags=5, selectlags="AIC")  )
summary( ur.df(dy,  type="trend", lags=5, selectlags="AIC")  )

## Testing for GARCH 
plot(dy2);      Sys.slepp(5)
acf(dy2);                   # (dy)^2 autocorrelated ?
LjungBox(dy2,10)

######################## y
### GARCH for y   -  Version 1
mod.arima <- arima(y, order=c(1,0,1), include.mean=TRUE)
coef(mod.arima);    mod.arima$aic; 
res <- mod.arima$residuals; 

plot(res);      Sys.sleep(5)  
acf(res);       Sys.sleep(5)  
plot(res^2);    Sys.sleep(5)
plot(abs(res)); Sys.sleep(5)
acf(res^2)
LjungBox(res^2,10)

# GARCH 
mm <- garch(res, order=c(0,1) )     # ARCH(1)
coef(mm);    aic <- 2*mm$n.likeli + 2*sum(mm$order); aic   # AIC --> min
#summary( garch(res, order=c(0,1) ) ) 

mm <- garch(res, order=c(1,1) )    # GARCH(1,1)
coef(mm);    aic <- 2*mm$n.likeli + 2*sum(mm$order); aic   # AIC
#summary( garch(res, order=c(1,1) ) ) 

#garch(res, order=c(1,0) )         # GARCH(1,0)  # refuses estimation

##  not stable
spec1 = ugarchspec(                      # here ARCH(1) == c(1,0) !!!
   variance.model = list( model = "sGARCH", garchOrder = c(1, 1) ),  
   mean.model = list( armaOrder = c(0, 0), include.mean = FALSE ) 
                  )               
fit = ugarchfit(spec=spec1,data=res)
#fit
coef(fit)             # head(sigma(fit))


### GARCH for y   -  Version 2   ## not stable
spec1 = ugarchspec(
   variance.model = list( model = "sGARCH", garchOrder = c(1, 0) ), 
   mean.model = list( armaOrder = c(1, 1), include.mean = TRUE ) 
                  )           
fit = ugarchfit(spec=spec1,data=y)    # fit
coef(fit)                             # head(sigma(fit))


################## dy

### GARCH for dy   -  Version 1
mod.arima <- arima(dy, order=c(0,0,1), include.mean=TRUE)
coef(mod.arima);    mod.arima$aic; 
res <- mod.arima$residuals; 

plot(res);      Sys.sleep(5)  
acf(res);       Sys.sleep(5)  
plot(res^2);    Sys.sleep(5)
plot(abs(res)); Sys.sleep(5)
acf(res^2)
LjungBox(res^2,10)

# GARCH 
mm <- garch(res, order=c(0,1) )     # ARCH(1)
coef(mm);    aic <- 2*mm$n.likeli + 2*sum(mm$order); aic   # AIC --> min
#summary( garch(res, order=c(0,1) ) ) 

mm <- garch(res, order=c(1,1) )    # GARCH(1,1)
coef(mm);    aic <- 2*mm$n.likeli + 2*sum(mm$order); aic   # AIC
#summary( garch(res, order=c(1,1) ) ) 


##
spec1 = ugarchspec(             # not useful
   variance.model = list( model = "sGARCH", garchOrder = c(0, 1) ),  
   mean.model = list( armaOrder = c(0, 0), include.mean = FALSE ) 
                  )               
fit = ugarchfit(spec=spec1,data=res)
#fit
coef(fit)             # head(sigma(fit))



### GARCH for dy   -  Version 2
spec1 = ugarchspec(
   variance.model = list( model = "sGARCH", garchOrder = c(1, 1) ), 
   mean.model = list( armaOrder = c(0, 1), include.mean = TRUE ) 
                  )
                 
fit = ugarchfit(spec=spec1,data=dy)
#fit
coef(fit)             # head(sigma(fit))










###################

spec1 = ugarchspec(
   variance.model = list( model = "sGARCH", garchOrder = c(1, 0) ), 
   mean.model = list( armaOrder = c(1, 1), include.mean = TRUE ) 
                  )
                 
fit = ugarchfit(spec=spec1,data=y)
#fit
coef(fit)             # head(sigma(fit))








ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1), 
submodel = NULL, external.regressors = NULL, variance.targeting = FALSE), 
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE, archm = FALSE, 
archpow = 1, arfima = FALSE, external.regressors = NULL, archex = FALSE), 
distribution.model = "norm", start.pars = list(), fixed.pars = list(), ...)