#=======================================================

# clear Global Environment
rm(list = ls())

install.packages("Hmisc")
install.packages("fGarch")
install.packages("tseries")
install.packages("rugarch") 
install.packages('knitr')     
install.packages('lmtest')      
install.packages('zoo') 
install.packages('car')   

library(Hmisc)
library(car)     # gia  tipika sfalmata tou White
library(zoo)     #  gia time series functions  
library(knitr)   # gia kable()
library(lmtest)  # gia `coeftest()` and `bptest()`.
library(fGarch)  # gia na ektimoume arch/garch montela
library(rugarch) #gia GARCH models

#=======================================================
# Load Hedge fund data
mydata<- read.table("C:/Vrontos/mathimata/A-mathimata-2020-2021/BSc-Econometrics/Vrontos/Ergastirio/Example2/R/example2_data.txt")
dim(mydata)
y <- mydata$V1
x1 <- mydata$V2
x2 <- mydata$V3
x3 <- mydata$V4
x4 <- mydata$V5
x5 <- mydata$V6
x6 <- mydata$V7
x7 <- mydata$V8
x8 <- mydata$V9

#=======================================================

# Create a time series object and plot
y=ts(y, frequency=1, start = c(1989,4))
y
plot(y, type="l", col="blue")

# Explore linear relationships
# Correlation coefficients
cor(cbind(y,x1,x2,x3,x4,x5,x6,x7,x8))
# Correlation coefficients and p-values
rcorr(as.matrix(cbind(y,x1,x2,x3,x4,x5,x6,x7,x8))) 
# Scatterplot of all variables
pairs(cbind(y,x1,x2,x3,x4,x5,x6,x7,x8)) 

#=======================================================

# 1a.RUN MULTIPLE REGRESSION MODEL SELECTION TECHNIQUES
# Fit the full model
fitall <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)
summary(fitall)

# Backward Elimination method
stepBE<-step(fitall, 
             scope=list(lower = ~ 1, upper = ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8),
             direction="backward")
stepBE

# Forward Selection method
fitnull<-lm(y ~ 1)
stepFS<-step(fitnull, 
             scope=list(lower = ~ 1, upper = ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8),
             direction="forward")
stepFS

# Stepwise Selection method
stepSS<-step(fitnull, 
             scope=list(lower = ~ 1, upper = ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8),
             direction="both")
stepSS

stepBE
stepFS
stepSS

#=======================================================

# BUT !!!!  what tell us the characteristics of the data??

# Summary Statistics and plots
plot(y, type="l", main="monthly returns")
hist(y, main="histogram of returns")
qqnorm(y,main="Normal QQplot of y") # normal Q-Q plot
qqline(y) 
acf(y, 24, main="ACF of returns")   
pacf(y, 24, main="PACF of returns") 
acf(y^2,50, main="ACF of squared returns")
pacf(y^2, 50, main="PACF of squared returns")

Box.test(y,lag=12,type="Box-Pierce") 
Box.test(y^2,lag=2,type="Ljung") 

# Normality Test
jarque.bera.test(y)
shapiro.test(y) 

#=======================================================

# 1. RUN MULTIPLE REGRESSION MODEL
MRres <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8)
summary(MRres)

ehat <-resid(MRres)
plot(ehat, type='l')
abline(h=0, lty=2)


# AUTOCORRELATION TESTS
#======================
# Autocorrelation of the residuals
acf(MRres$residuals, 36)
pacf(MRres$residuals, 36)

# autocorrelation plot
correlogram <- acf(ehat)
table_correlog<-data.frame(correlogram$lag, correlogram$acf)
table_correlog

# Box-Pierce test
Box.test(ehat,12)
Box.test(ehat,12, type = "Ljung")
Q_stat<-Box.test(ehat, lag = 12, type = "Ljung")
Q_stat

for(k in 1:5){
  print(Box.test(ehat, lag = k, type = "Ljung-Box"))
}

# Durbin-Watson test
dwtest(MRres)  

# Breusch-Godfrey test
a <- bgtest(MRres, order=1, type="Chisq")
a
b <- bgtest(MRres, order=2, type="Chisq")
b
dfr <- data.frame(rbind(a[c(1,2,4)], b[c(1,2,4)]))
names(dfr)<-c(  "LM Statistic", "Lags", "p-Value")
kable(dfr, caption="Breusch-Godfrey test for our model")


# HETEROSKEDASTICITY TESTS
#=========================
ehat <- residuals(MRres)
ehat2 <- residuals(MRres)^2
ehat_abs <- abs(residuals(MRres))
yhat <- fitted(MRres)

plot(yhat,ehat, xlab="fitted values", ylab="residuals")
abline(h=0, col="blue")
plot(yhat,ehat2, xlab="fitted values", ylab="squared residuals")
plot(yhat,ehat_abs, xlab="fitted values", ylab="absolute residuals")

# Autocorrelation of the squared residuals
acf(MRres$residuals^2, 36)
pacf(MRres$residuals^2, 36)

#Breusch-Pagan-Godfrey test
# First way
bptest(MRres)

# Second way
alpha <- 0.05
# Auxilliary Regression:
aux_re <- lm(ehat2 ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8, data=mydata)
N <- nobs(aux_re)   # sample size
S <-aux_re$rank - 1 # number of explanatory variables (without constant)

summary_aux_re <- summary(aux_re)
R_squared <- summary_aux_re$r.squared  
LM_stat <- N*R_squared 
LM_stat

#Chi-square critical value
Chi_sq_crit <- qchisq(1-alpha, S)
Chi_sq_crit
p_value <- 1-pchisq(LM_stat ,S)
p_value


# White test
# First way
bptest(MRres, ~ x1+x2+x3+x4+x5+x6+x7+x8+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+I(x7^2)+I(x8^2))

# Second way
aux_re_2 <- lm(ehat2 ~ x1+x2+x3+x4+x5+x6+x7+x8+I(x1^2)+I(x2^2)+I(x3^2)+I(x4^2)+I(x5^2)+I(x6^2)+I(x7^2)+I(x8^2), data=mydata)
N <- nobs(aux_re_2)   # sample size
S <-aux_re_2$rank - 1 # number of explanatory variables (without constant)
summary_aux_re_2 <- summary(aux_re_2)
R_squared <- summary_aux_re_2$r.squared  
LM_stat <- N*R_squared
LM_stat
p_value <- 1-pchisq(LM_stat, S)
p_value


# NORMALITY TESTS
#================
jarque.bera.test(MRres$residuals)
shapiro.test(MRres$residuals) 
qqnorm(MRres$residuals) 
qqline(MRres$residuals)  


#=======================================================

# AR1 for regression residuals
ar1res=arima(residuals(MRres), order=c(1,0,0), include.mean=FALSE)
ar1res$coef
# Autocorrelation of the AR(1) residuals
acf(residuals(ar1res),36)
pacf(residuals(ar1res),36)
# Autocorrelation of the AR(1) squared residuals
acf(residuals(ar1res)^2, 36)
pacf(residuals(ar1res)^2, 36)
# Normality test of the AR(1)
jarque.bera.test(residuals(ar1res))
shapiro.test(residuals(ar1res)) 
qqnorm(residuals(ar1res)) 
qqline(residuals(ar1res))  

#=======================================================

# 2. RUN MULTIPLE REGRESSION MODEL with time series components
MRAR1res = arima(y, order=c(1,0,0), xreg=cbind(x1,x2,x3,x4,x5,x6,x7,x8)) 
MRAR1res
# Diagnostic tests for the residuals
# Autocorrelation of the residuals
acf(MRAR1res$residuals, 36)
pacf(MRAR1res$residuals, 36)
# Autocorrelation of the squared residuals
acf(MRAR1res$residuals^2, 36)
pacf(MRAR1res$residuals^2, 36)
# Normality test
jarque.bera.test(MRAR1res$residuals)
shapiro.test(MRAR1res$residuals) 
qqnorm(MRAR1res$residuals) 
qqline(MRAR1res$residuals)  

#=======================================================

# 3. RUN MULTIPLE REGRESSION MODEL + ARMA + GARCH
X<-matrix(cbind(x1,x2,x3,x4,x5,x6,x7,x8),ncol=8)
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder=c(1,1)), 
                   mean.model = list(armaOrder=c(1,0), include.mean = TRUE, external.regressors = X),
                   distribution.model = "norm")
spec
modelres <- ugarchfit(spec = spec, data = y)
modelres

pred <- ugarchforecast(modelres, 
                       n.ahead = 6, 
                       external.forecasts =list(mean(x1),mean(x2),
                                                mean(x3),mean(x4),
                                                mean(x5),mean(x6),
                                                mean(x7),mean(x8)))
pred

#=======================================================