library(PerformanceAnalytics)
library(quantmod)
library(rugarch)
library(car)
library(FinTS)
library("rmutil")
library("normtest")
library("tseries")
library(urca)
options(digits=4)
require(ggplot2)


################################################################################################
#####                                Simulacni studie

###############  GARCH(1,1) Normal 100
###############################################################################################

spec1norm <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                        fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                          "beta1" = 0.85))

show(spec1norm)

pocetSim=10000

j=100

zamitSBT = 0
zamitNSBT = 0
zamitPSBT = 0
zamitJoint = 0

omega100 <- c()
alpha1100 <- c()
beta1100 <- c()

kounter = 0
i=1

spec_norm <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE))  


repeat{
  
  x_obj_norm <- ugarchpath(spec1norm, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm <- x_obj_norm@path$seriesSim
  
  try({
    
    fit_norm <- ugarchfit(spec_norm, data = x1norm)
    
    omega100[i] <- coef(fit_norm)[1]
    alpha1100[i] <- coef(fit_norm)[2]
    beta1100[i] <- coef(fit_norm)[3]
    
    
    #####Sign bias tests
    
    if(is(fit_norm, "uGARCHfilter")) z = fit_norm@filter$z else z = z = fit_norm@fit$z
    res = as.numeric(residuals(fit_norm))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_norm)

    SBTp <- summary(fitSB)$coefficients[2,4]
    NSBTp <- summary(fitNSB)$coefficients[2,4]
    PSBTp <- summary(fitPSB)$coefficients[2,4]
    JointTp <- signb$prob[4]
    
    if (SBTp<0.05) {zamitSBT=zamitSBT+1}
    if (NSBTp<0.05) {zamitNSBT=zamitNSBT+1}
    if (PSBTp<0.05) {zamitPSBT=zamitPSBT+1}
    if (JointTp<0.05) {zamitJoint=zamitJoint+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT=zamitSBT/pocetSim*100
ratioNSBT=zamitNSBT/pocetSim*100
ratioPSBT=zamitPSBT/pocetSim*100
ratioJoint=zamitJoint/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio <- c(ratioSBT,ratioNSBT,ratioPSBT,ratioJoint)

ratiomatice <- rbind(testy,ratio)

print(j)
print(ratiomatice)

print(kounter)

print("omega")
print(summary(omega100))
print(sd(omega100))

print("alpha1")
print(summary(alpha1100))
print(sd(alpha1100))

print("beta1")
print(summary(beta1100))
print(sd(beta1100))
###############################################################################################

###############  GARCH(1,1) Normal 500
###############################################################################################

spec1norm <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                        fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                          "beta1" = 0.85))

show(spec1norm)

pocetSim=10000

j=500

zamitSBT = 0
zamitNSBT = 0
zamitPSBT = 0
zamitJoint = 0

omega500 <- c()
alpha1500 <- c()
beta1500 <- c()

kounter = 0
i=1

spec_norm <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE))  


repeat{
  
  x_obj_norm <- ugarchpath(spec1norm, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm <- x_obj_norm@path$seriesSim
  
  try({
    
    fit_norm <- ugarchfit(spec_norm, data = x1norm)
    
    omega500[i] <- coef(fit_norm)[1]
    alpha1500[i] <- coef(fit_norm)[2]
    beta1500[i] <- coef(fit_norm)[3]
    
    
    #####Sign bias tests
    
    if(is(fit_norm, "uGARCHfilter")) z = fit_norm@filter$z else z = z = fit_norm@fit$z
    res = as.numeric(residuals(fit_norm))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_norm)

    SBTp <- summary(fitSB)$coefficients[2,4]
    NSBTp <- summary(fitNSB)$coefficients[2,4]
    PSBTp <- summary(fitPSB)$coefficients[2,4]
    JointTp <- signb$prob[4]
    
    if (SBTp<0.05) {zamitSBT=zamitSBT+1}
    if (NSBTp<0.05) {zamitNSBT=zamitNSBT+1}
    if (PSBTp<0.05) {zamitPSBT=zamitPSBT+1}
    if (JointTp<0.05) {zamitJoint=zamitJoint+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT=zamitSBT/pocetSim*100
ratioNSBT=zamitNSBT/pocetSim*100
ratioPSBT=zamitPSBT/pocetSim*100
ratioJoint=zamitJoint/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio <- c(ratioSBT,ratioNSBT,ratioPSBT,ratioJoint)

ratiomatice <- rbind(testy,ratio)

print(j)
print(ratiomatice)

print(kounter)

print("omega")
print(summary(omega500))
print(sd(omega500))

print("alpha1")
print(summary(alpha1500))
print(sd(alpha1500))

print("beta1")
print(summary(beta1500))
print(sd(beta1500))
###############################################################################################

###############  GARCH(1,1) Normal 1000
###############################################################################################

spec1norm <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                        fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                          "beta1" = 0.85))

show(spec1norm)

pocetSim=10000

j=1000

zamitSBT = 0
zamitNSBT = 0
zamitPSBT = 0
zamitJoint = 0

omega1000 <- c()
alpha11000 <- c()
beta11000 <- c()

kounter = 0
i=1

spec_norm <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE))  


repeat{
  
  x_obj_norm <- ugarchpath(spec1norm, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm <- x_obj_norm@path$seriesSim
  
  try({
    
    fit_norm <- ugarchfit(spec_norm, data = x1norm)
    
    omega1000[i] <- coef(fit_norm)[1]
    alpha11000[i] <- coef(fit_norm)[2]
    beta11000[i] <- coef(fit_norm)[3]
    
    
    #####Sign bias tests
    
    if(is(fit_norm, "uGARCHfilter")) z = fit_norm@filter$z else z = z = fit_norm@fit$z
    res = as.numeric(residuals(fit_norm))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_norm)

    SBTp <- summary(fitSB)$coefficients[2,4]
    NSBTp <- summary(fitNSB)$coefficients[2,4]
    PSBTp <- summary(fitPSB)$coefficients[2,4]
    JointTp <- signb$prob[4]
    
    if (SBTp<0.05) {zamitSBT=zamitSBT+1}
    if (NSBTp<0.05) {zamitNSBT=zamitNSBT+1}
    if (PSBTp<0.05) {zamitPSBT=zamitPSBT+1}
    if (JointTp<0.05) {zamitJoint=zamitJoint+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT=zamitSBT/pocetSim*100
ratioNSBT=zamitNSBT/pocetSim*100
ratioPSBT=zamitPSBT/pocetSim*100
ratioJoint=zamitJoint/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio <- c(ratioSBT,ratioNSBT,ratioPSBT,ratioJoint)

ratiomatice <- rbind(testy,ratio)

print(j)
print(ratiomatice)

print(kounter)

print("omega")
print(summary(omega1000))
print(sd(omega1000))

print("alpha1")
print(summary(alpha11000))
print(sd(alpha11000))

print("beta1")
print(summary(beta11000))
print(sd(beta11000))
###############################################################################################

###############  GARCH(1,1) Laplace 100 
###############################################################################################

spec1_lapl <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                         mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                         distribution.model = "ged",
                         fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                           "beta1" = 0.85, "shape"=1))

show(spec1_lapl)

pocetSim=10000

j=100

zamitSBT_lapl = 0
zamitNSBT_lapl = 0
zamitPSBT_lapl = 0
zamitJoint_lapl = 0

omega100_l <- c()
alpha1100_l <- c()
beta1100_l <- c()

kounter = 0
i=1

spec_lapl <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        distribution.model = "ged",
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                        fixed.pars = list("shape"=1))

repeat{
  
  x_obj_lapl <- ugarchpath(spec1_lapl, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl <- x_obj_lapl@path$seriesSim
  
  try({
    
    fit_lapl <- ugarchfit(spec_lapl, data = x1_lapl)
    
    omega100_l[i] <- coef(fit_lapl)[1]
    alpha1100_l[i] <- coef(fit_lapl)[2]
    beta1100_l[i] <- coef(fit_lapl)[3]

    #####Sign bias tests
    
    if(is(fit_lapl, "uGARCHfilter")) z = fit_lapl@filter$z else z = z = fit_lapl@fit$z
    res = as.numeric(residuals(fit_lapl))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)
  
    signb <- signbias(fit_lapl)

    SBTp_lapl <- summary(fitSB)$coefficients[2,4]
    NSBTp_lapl <- summary(fitNSB)$coefficients[2,4]
    PSBTp_lapl <- summary(fitPSB)$coefficients[2,4]
    JointTp_lapl <- signb$prob[4]
    
    if (SBTp_lapl<0.05) {zamitSBT_lapl=zamitSBT_lapl+1}
    if (NSBTp_lapl<0.05) {zamitNSBT_lapl=zamitNSBT_lapl+1}
    if (PSBTp_lapl<0.05) {zamitPSBT_lapl=zamitPSBT_lapl+1}
    if (JointTp_lapl<0.05) {zamitJoint_lapl=zamitJoint_lapl+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_lapl=zamitSBT_lapl/pocetSim*100
ratioNSBT_lapl=zamitNSBT_lapl/pocetSim*100
ratioPSBT_lapl=zamitPSBT_lapl/pocetSim*100
ratioJoint_lapl=zamitJoint_lapl/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_lapl <- c(ratioSBT_lapl,ratioNSBT_lapl,ratioPSBT_lapl,ratioJoint_lapl)

ratiomatice_lapl <- rbind(testy,ratio_lapl)

print(j)
print(ratiomatice_lapl)

print(kounter)

print("omega")
print(summary(omega100_l))
print(sd(omega100))

print("alpha1")
print(summary(alpha1100_l))
print(sd(alpha1100))

print("beta1")
print(summary(beta1100_l))
print(sd(beta1100))
###############################################################################################

###############  GARCH(1,1) Laplace 500
###############################################################################################

spec1_lapl <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                         mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                         distribution.model = "ged",
                         fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                           "beta1" = 0.85, "shape"=1))

show(spec1_lapl)

pocetSim=10000

j=500

zamitSBT_lapl = 0
zamitNSBT_lapl = 0
zamitPSBT_lapl = 0
zamitJoint_lapl = 0

omega500_l <- c()
alpha1500_l <- c()
beta1500_l <- c()

kounter = 0
i=1

spec_lapl <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        distribution.model = "ged",
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                        fixed.pars = list("shape"=1))

repeat{
  
  x_obj_lapl <- ugarchpath(spec1_lapl, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl <- x_obj_lapl@path$seriesSim
  
  try({
    
    fit_lapl <- ugarchfit(spec_lapl, data = x1_lapl)
    
    omega500_l[i] <- coef(fit_lapl)[1]
    alpha1500_l[i] <- coef(fit_lapl)[2]
    beta1500_l[i] <- coef(fit_lapl)[3]

    #####Sign bias tests
    
    if(is(fit_lapl, "uGARCHfilter")) z = fit_lapl@filter$z else z = z = fit_lapl@fit$z
    res = as.numeric(residuals(fit_lapl))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_lapl)

    SBTp_lapl <- summary(fitSB)$coefficients[2,4]
    NSBTp_lapl <- summary(fitNSB)$coefficients[2,4]
    PSBTp_lapl <- summary(fitPSB)$coefficients[2,4]
    JointTp_lapl <- signb$prob[4]
    
    if (SBTp_lapl<0.05) {zamitSBT_lapl=zamitSBT_lapl+1}
    if (NSBTp_lapl<0.05) {zamitNSBT_lapl=zamitNSBT_lapl+1}
    if (PSBTp_lapl<0.05) {zamitPSBT_lapl=zamitPSBT_lapl+1}
    if (JointTp_lapl<0.05) {zamitJoint_lapl=zamitJoint_lapl+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_lapl=zamitSBT_lapl/pocetSim*100
ratioNSBT_lapl=zamitNSBT_lapl/pocetSim*100
ratioPSBT_lapl=zamitPSBT_lapl/pocetSim*100
ratioJoint_lapl=zamitJoint_lapl/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_lapl <- c(ratioSBT_lapl,ratioNSBT_lapl,ratioPSBT_lapl,ratioJoint_lapl)

ratiomatice_lapl <- rbind(testy,ratio_lapl)

print(j)
print(ratiomatice_lapl)

print(kounter)

print("omega")
print(summary(omega500_l))
print(sd(omega500))

print("alpha1")
print(summary(alpha1500_l))
print(sd(alpha1500))

print("beta1")
print(summary(beta1500_l))
print(sd(beta1500))
###############################################################################################

###############  GARCH(1,1) Laplace 1000
###############################################################################################

spec1_lapl <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                         mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                         distribution.model = "ged",
                         fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                           "beta1" = 0.85, "shape"=1))

show(spec1_lapl)

pocetSim=10000

j=1000

zamitSBT_lapl = 0
zamitNSBT_lapl = 0
zamitPSBT_lapl = 0
zamitJoint_lapl = 0

omega1000_l <- c()
alpha11000_l <- c()
beta11000_l <- c()

kounter = 0
i=1

spec_lapl <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                        distribution.model = "ged",
                        mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                        fixed.pars = list("shape"=1))

repeat{
  
  x_obj_lapl <- ugarchpath(spec1_lapl, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl <- x_obj_lapl@path$seriesSim
  
  try({
    
    fit_lapl <- ugarchfit(spec_lapl, data = x1_lapl)
    
    omega1000_l[i] <- coef(fit_lapl)[1]
    alpha11000_l[i] <- coef(fit_lapl)[2]
    beta11000_l[i] <- coef(fit_lapl)[3]

    #####Sign bias tests
    
    if(is(fit_lapl, "uGARCHfilter")) z = fit_lapl@filter$z else z = z = fit_lapl@fit$z
    res = as.numeric(residuals(fit_lapl))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_lapl)

    SBTp_lapl <- summary(fitSB)$coefficients[2,4]
    NSBTp_lapl <- summary(fitNSB)$coefficients[2,4]
    PSBTp_lapl <- summary(fitPSB)$coefficients[2,4]
    JointTp_lapl <- signb$prob[4]
    
    if (SBTp_lapl<0.05) {zamitSBT_lapl=zamitSBT_lapl+1}
    if (NSBTp_lapl<0.05) {zamitNSBT_lapl=zamitNSBT_lapl+1}
    if (PSBTp_lapl<0.05) {zamitPSBT_lapl=zamitPSBT_lapl+1}
    if (JointTp_lapl<0.05) {zamitJoint_lapl=zamitJoint_lapl+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_lapl=zamitSBT_lapl/pocetSim*100
ratioNSBT_lapl=zamitNSBT_lapl/pocetSim*100
ratioPSBT_lapl=zamitPSBT_lapl/pocetSim*100
ratioJoint_lapl=zamitJoint_lapl/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_lapl <- c(ratioSBT_lapl,ratioNSBT_lapl,ratioPSBT_lapl,ratioJoint_lapl)

ratiomatice_lapl <- rbind(testy,ratio_lapl)

print(j)
print(ratiomatice_lapl)

print(kounter)

print("omega")
print(summary(omega1000_l))
print(sd(omega1000))

print("alpha1")
print(summary(alpha11000_l))
print(sd(alpha11000))

print("beta1")
print(summary(beta11000_l))
print(sd(beta11000))
###############################################################################################

###############  GARCH(1,1) t 100
###############################################################################################

spec1t <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                     mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                     distribution.model = "std",
                     fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                       "beta1" = 0.85, "shape"=4))

show(spec1t)

pocetSim=10000

j=100

zamitSBT_t = 0
zamitNSBT_t = 0
zamitPSBT_t = 0
zamitJoint_t = 0

omega100_t <- c()
alpha1100_t <- c()
beta1100_t <- c()

kounter=0
i=1

spec_t <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                     distribution.model = "std",
                     mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                     fixed.pars = list("shape"=4))



repeat{
  x_obj_t <- ugarchpath(spec1t, n.sim = j, n.start = 1000, m.sim=1)
  
  x1t <- x_obj_t@path$seriesSim
  
  try({
    
    fit_t <- ugarchfit(spec_t, data = x1t)
    
    omega100_t[i] <- coef(fit_t)[1]
    alpha1100_t[i] <- coef(fit_t)[2]
    beta1100_t[i] <- coef(fit_t)[3]
    
    #####Sign bias tests
    
    if(is(fit_t, "uGARCHfilter")) z = fit_t@filter$z else z = z = fit_t@fit$z
    res = as.numeric(residuals(fit_t))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_t)

    SBTp_t <- summary(fitSB)$coefficients[2,4]
    NSBTp_t <- summary(fitNSB)$coefficients[2,4]
    PSBTp_t <- summary(fitPSB)$coefficients[2,4]
    JointTp_t <- signb$prob[4]
    
    if (SBTp_t<0.05) {zamitSBT_t=zamitSBT_t+1}
    if (NSBTp_t<0.05) {zamitNSBT_t=zamitNSBT_t+1}
    if (PSBTp_t<0.05) {zamitPSBT_t=zamitPSBT_t+1}
    if (JointTp_t<0.05) {zamitJoint_t=zamitJoint_t+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_t=zamitSBT_t/pocetSim*100
ratioNSBT_t=zamitNSBT_t/pocetSim*100
ratioPSBT_t=zamitPSBT_t/pocetSim*100
ratioJoint_t=zamitJoint_t/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_t <- c(ratioSBT_t,ratioNSBT_t,ratioPSBT_t,ratioJoint_t)

ratiomatice_t <- rbind(testy,ratio_t)

print(j)
print(ratiomatice_t)

print(kounter)

print("omega")
print(summary(omega100_t))
print(sd(omega100_t))

print("alpha1")
print(summary(alpha1100_t))
print(sd(alpha1100_t))

print("beta1")
print(summary(beta1100_t))
print(sd(beta1100_t))
###############################################################################################

###############  GARCH(1,1) t 500
###############################################################################################

spec1t <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                     mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                     distribution.model = "std",
                     fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                       "beta1" = 0.85, "shape"=4))

show(spec1t)

pocetSim=10000

j=500

zamitSBT_t = 0
zamitNSBT_t = 0
zamitPSBT_t = 0
zamitJoint_t = 0

omega500_t <- c()
alpha1500_t <- c()
beta1500_t <- c()

kounter=0
i=1

spec_t <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                     distribution.model = "std",
                     mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                     fixed.pars = list("shape"=4))



repeat{
  x_obj_t <- ugarchpath(spec1t, n.sim = j, n.start = 1000, m.sim=1)
  
  x1t <- x_obj_t@path$seriesSim
  
  try({
    
    fit_t <- ugarchfit(spec_t, data = x1t)
    
    omega500_t[i] <- coef(fit_t)[1]
    alpha1500_t[i] <- coef(fit_t)[2]
    beta1500_t[i] <- coef(fit_t)[3]
    
    #####Sign bias tests
    
    if(is(fit_t, "uGARCHfilter")) z = fit_t@filter$z else z = z = fit_t@fit$z
    res = as.numeric(residuals(fit_t))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_t)

    SBTp_t <- summary(fitSB)$coefficients[2,4]
    NSBTp_t <- summary(fitNSB)$coefficients[2,4]
    PSBTp_t <- summary(fitPSB)$coefficients[2,4]
    JointTp_t <- signb$prob[4]
    
    if (SBTp_t<0.05) {zamitSBT_t=zamitSBT_t+1}
    if (NSBTp_t<0.05) {zamitNSBT_t=zamitNSBT_t+1}
    if (PSBTp_t<0.05) {zamitPSBT_t=zamitPSBT_t+1}
    if (JointTp_t<0.05) {zamitJoint_t=zamitJoint_t+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_t=zamitSBT_t/pocetSim*100
ratioNSBT_t=zamitNSBT_t/pocetSim*100
ratioPSBT_t=zamitPSBT_t/pocetSim*100
ratioJoint_t=zamitJoint_t/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_t <- c(ratioSBT_t,ratioNSBT_t,ratioPSBT_t,ratioJoint_t)

ratiomatice_t <- rbind(testy,ratio_t)

print(j)
print(ratiomatice_t)

print(kounter)

print("omega")
print(summary(omega500_t))
print(sd(omega500_t))

print("alpha1")
print(summary(alpha1500_t))
print(sd(alpha1500_t))

print("beta1")
print(summary(beta1500_t))
print(sd(beta1500_t))
###############################################################################################

###############  GARCH(1,1) t 1000
###############################################################################################

spec1t <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                     mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                     distribution.model = "std",
                     fixed.pars = list("omega" = 0.05, "alpha1" = 0.1,
                                       "beta1" = 0.85, "shape"=4))

show(spec1t)

pocetSim=10000

j=1000

zamitSBT_t = 0
zamitNSBT_t = 0
zamitPSBT_t = 0
zamitJoint_t = 0

omega1000_t <- c()
alpha11000_t <- c()
beta11000_t <- c()

kounter=0
i=1

spec_t <- ugarchspec(variance.model = list(garchOrder=c(1,1)),
                     distribution.model = "std",
                     mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                     fixed.pars = list("shape"=4))



repeat{
  x_obj_t <- ugarchpath(spec1t, n.sim = j, n.start = 1000, m.sim=1)
  
  x1t <- x_obj_t@path$seriesSim
  
  try({
    
    fit_t <- ugarchfit(spec_t, data = x1t)
    
    omega1000_t[i] <- coef(fit_t)[1]
    alpha11000_t[i] <- coef(fit_t)[2]
    beta11000_t[i] <- coef(fit_t)[3]
    
    #####Sign bias tests
    
    if(is(fit_t, "uGARCHfilter")) z = fit_t@filter$z else z = z = fit_t@fit$z
    res = as.numeric(residuals(fit_t))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])

    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_t)

    SBTp_t <- summary(fitSB)$coefficients[2,4]
    NSBTp_t <- summary(fitNSB)$coefficients[2,4]
    PSBTp_t <- summary(fitPSB)$coefficients[2,4]
    JointTp_t <- signb$prob[4]
    
    if (SBTp_t<0.05) {zamitSBT_t=zamitSBT_t+1}
    if (NSBTp_t<0.05) {zamitNSBT_t=zamitNSBT_t+1}
    if (PSBTp_t<0.05) {zamitPSBT_t=zamitPSBT_t+1}
    if (JointTp_t<0.05) {zamitJoint_t=zamitJoint_t+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_t=zamitSBT_t/pocetSim*100
ratioNSBT_t=zamitNSBT_t/pocetSim*100
ratioPSBT_t=zamitPSBT_t/pocetSim*100
ratioJoint_t=zamitJoint_t/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_t <- c(ratioSBT_t,ratioNSBT_t,ratioPSBT_t,ratioJoint_t)

ratiomatice_t <- rbind(testy,ratio_t)

print(j)
print(ratiomatice_t)

print(kounter)

print("omega")
print(summary(omega1000_t))
print(sd(omega1000_t))

print("alpha1")
print(summary(alpha11000_t))
print(sd(alpha11000_t))

print("beta1")
print(summary(beta11000_t))
print(sd(beta11000_t))
###############################################################################################

###############  EGARCH(1,1) Normal 100 odhad parametru
###############################################################################################
spec1norm_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                               fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                 "beta1" = 0.7,"gamma1"=0.28))

show(spec1norm_egarch)

pocetSim=10000

j=100

omega100_egarch <- c()
alpha1100_egarch <- c()
beta1100_egarch <- c()
gamma100_egarch <-c()

kounter=0
i=1

spec_norm_egarch <- ugarchspec(variance.model = list(model="eGARCH",
                    garchOrder=c(1,1)),
                    mean.model = list(armaOrder = c(0,0), include.mean = FALSE))

repeat {
  
  x_obj_norm_egarch <- ugarchpath(spec1norm_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm_egarch <- x_obj_norm_egarch@path$seriesSim
  
  
  try({
    
    fit_norm_egarch <- ugarchfit(spec_norm_egarch, data = x1norm_egarch)
    
    omega100_egarch[i] <- coef(fit_norm_egarch)[1]
    alpha1100_egarch[i] <- coef(fit_norm_egarch)[2]
    beta1100_egarch[i] <- coef(fit_norm_egarch)[3]
    gamma100_egarch[i] <- coef(fit_norm_egarch)[4]
    
    kounter=kounter+1
    
    i=i+1
    
    if (i%%1000==0){print(i)}
    
  } , silent = TRUE)
  
  if (i==pocetSim+1) {break}
}
print(j)

print(kounter)

print("omega")
omega100_egarch<-omega100_egarch[!is.na(omega100_egarch)]
print(summary(omega100_egarch))
print(sd(omega100_egarch))

print("alpha1")
alpha1100_egarch<-alpha1100_egarch[!is.na(alpha1100_egarch)]
print(summary(alpha1100_egarch))
print(sd(alpha1100_egarch))

print("beta1")
beta1100_egarch<-beta1100_egarch[!is.na(beta1100_egarch)]
print(summary(beta1100_egarch))
print(sd(beta1100_egarch))

print("gamma")
gamma100_egarch<-gamma100_egarch[!is.na(gamma100_egarch)]
print(summary(gamma100_egarch))
print(sd(gamma100_egarch))
###############################################################################################

###############  EGARCH(1,1) Normal 500 odhad parametru
###############################################################################################
spec1norm_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                               fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                 "beta1" = 0.7,"gamma1"=0.28))

show(spec1norm_egarch)

pocetSim=10000

j=500

omega500_egarch <- c()
alpha1500_egarch <- c()
beta1500_egarch <- c()
gamma500_egarch <-c()

kounter=0
i=1

spec_norm_egarch <- ugarchspec(variance.model = list(model="eGARCH",
                    garchOrder=c(1,1)),
                    mean.model = list(armaOrder = c(0,0), include.mean = FALSE))

repeat {
  
  x_obj_norm_egarch <- ugarchpath(spec1norm_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm_egarch <- x_obj_norm_egarch@path$seriesSim
  
  
  try({
    
    fit_norm_egarch <- ugarchfit(spec_norm_egarch, data = x1norm_egarch)
    
    omega500_egarch[i] <- coef(fit_norm_egarch)[1]
    alpha1500_egarch[i] <- coef(fit_norm_egarch)[2]
    beta1500_egarch[i] <- coef(fit_norm_egarch)[3]
    gamma500_egarch[i] <- coef(fit_norm_egarch)[4]
    
    kounter=kounter+1

    i=i+1

  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}

}
print(j)

print(kounter)

print("omega")
omega500_egarch<-omega500_egarch[!is.na(omega500_egarch)]
print(summary(omega500_egarch))
print(sd(omega500_egarch))

print("alpha1")
alpha1500_egarch<-alpha1500_egarch[!is.na(alpha1500_egarch)]
print(summary(alpha1500_egarch))
print(sd(alpha1500_egarch))

print("beta1")
beta1500_egarch<-beta1500_egarch[!is.na(beta1500_egarch)]
print(summary(beta1500_egarch))
print(sd(beta1500_egarch))

print("gamma")
gamma500_egarch<-gamma500_egarch[!is.na(gamma500_egarch)]
print(summary(gamma500_egarch))
print(sd(gamma500_egarch))
###############################################################################################

###############  EGARCH(1,1) Normal 1000 odhad parametru
###############################################################################################
spec1norm_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                               fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                 "beta1" = 0.7,"gamma1"=0.28))

show(spec1norm_egarch)

pocetSim=10000

j=1000

kounter=0
i=1

spec_norm_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE))

repeat {
  
  x_obj_norm_egarch <- ugarchpath(spec1norm_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm_egarch <- x_obj_norm_egarch@path$seriesSim
  
  
  try({
    
    fit_norm_egarch <- ugarchfit(spec_norm_egarch, data = x1norm_egarch)
    
    omega1000_egarch[i] <- coef(fit_norm_egarch)[1]
    alpha11000_egarch[i] <- coef(fit_norm_egarch)[2]
    beta11000_egarch[i] <- coef(fit_norm_egarch)[3]
    gamma1000_egarch[i] <- coef(fit_norm_egarch)[4]
    
    kounter=kounter+1

    i=i+1

  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}

}
print(j)

print(kounter)

print("omega")
omega1000_egarch<-omega1000_egarch[!is.na(omega1000_egarch)]
print(summary(omega1000_egarch))
print(sd(omega1000_egarch))

print("alpha1")
alpha11000_egarch<-alpha11000_egarch[!is.na(alpha11000_egarch)]
print(summary(alpha11000_egarch))
print(sd(alpha11000_egarch))

print("beta1")
beta11000_egarch<-beta11000_egarch[!is.na(beta11000_egarch)]
print(summary(beta11000_egarch))
print(sd(beta11000_egarch))

print("gamma")
gamma1000_egarch<-gamma1000_egarch[!is.na(gamma1000_egarch)]
print(summary(gamma1000_egarch))
print(sd(gamma1000_egarch))
###############################################################################################

###############  EGARCH(1,1) Laplace 100 odhad parametru
###############################################################################################
spec1_lapl_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                    garchOrder=c(1,1)),
                                mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                                distribution.model = "ged",
                                fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                  "beta1" = 0.7,"gamma1"=0.28, "shape"=1))

show(spec1_lapl_egarch)

pocetSim=10000

j=100

kounter=0  
i=1

spec_lapl_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "ged",
  fixed.pars = list("shape"=1))

repeat{
  
  x_obj_lapl_egarch <- ugarchpath(spec1_lapl_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl_egarch <- x_obj_lapl_egarch@path$seriesSim
  
  try({
    
    fit_lapl_egarch <- ugarchfit(spec_lapl_egarch, data = x1_lapl_egarch)
    
    omega100_egarch_l[i] <- coef(fit_lapl_egarch)[1]
    alpha1100_egarch_l[i] <- coef(fit_lapl_egarch)[2]
    beta1100_egarch_l[i] <- coef(fit_lapl_egarch)[3]
    gamma100_egarch_l[i] <- coef(fit_lapl_egarch)[4]
    
    kounter=kounter+1
    i=i+1

  } , silent = TRUE)

  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
  
}
print(j)

print(kounter)

print("omega")
print(summary(omega100_egarch_l))
print(sd(omega100_egarch_l))

print("alpha1")
print(summary(alpha1100_egarch_l))
print(sd(alpha1100_egarch_l))

print("beta1")
print(summary(beta1100_egarch_l))
print(sd(beta1100_egarch_l))

print("gamma")
print(summary(gamma100_egarch_l))
print(sd(gamma100_egarch_l))
###############################################################################################

###############  EGARCH(1,1) Laplace 500 odhad parametru
###############################################################################################
spec1_lapl_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                    garchOrder=c(1,1)),
                                mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                                distribution.model = "ged",
                                fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                  "beta1" = 0.7,"gamma1"=0.28, "shape"=1))

show(spec1_lapl_egarch)

pocetSim=10000

j=500

kounter=0  
i=1

spec_lapl_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "ged",
  fixed.pars = list("shape"=1))


repeat{
  x_obj_lapl_egarch <- ugarchpath(spec1_lapl_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl_egarch <- x_obj_lapl_egarch@path$seriesSim
  
  try({
    
    fit_lapl_egarch <- ugarchfit(spec_lapl_egarch, data = x1_lapl_egarch)
    
    omega500_egarch_l[i] <- coef(fit_lapl_egarch)[1]
    alpha1500_egarch_l[i] <- coef(fit_lapl_egarch)[2]
    beta1500_egarch_l[i] <- coef(fit_lapl_egarch)[3]
    gamma500_egarch_l[i] <- coef(fit_lapl_egarch)[4]
    
    kounter=kounter+1
    i=i+1
    
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
  
}
print(j)

print(kounter)

print("omega")
print(summary(omega500_egarch_l))
print(sd(omega500_egarch_l))

print("alpha1")
print(summary(alpha1500_egarch_l))
print(sd(alpha1500_egarch_l))

print("beta1")
print(summary(beta1500_egarch_l))
print(sd(beta1500_egarch_l))

print("gamma")
print(summary(gamma500_egarch_l))
print(sd(gamma500_egarch_l))
###############################################################################################

###############  EGARCH(1,1) Laplace 1000 odhad parametru
###############################################################################################
spec1_lapl_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                    garchOrder=c(1,1)),
                                mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                                distribution.model = "ged",
                                fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                  "beta1" = 0.7,"gamma1"=0.28, "shape"=1))

show(spec1_lapl_egarch)

pocetSim=10000

j=1000

omega1000_egarch_l <- c()
alpha11000_egarch_l <- c()
beta11000_egarch_l <- c()
gamma1000_egarch_l<-c()

kounter=0  
i=1

spec_lapl_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "ged",
  fixed.pars = list("shape"=1))


repeat{
  x_obj_lapl_egarch <- ugarchpath(spec1_lapl_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl_egarch <- x_obj_lapl_egarch@path$seriesSim
  
  try({
    
    fit_lapl_egarch <- ugarchfit(spec_lapl_egarch, data = x1_lapl_egarch)
    
    omega1000_egarch_l[i] <- coef(fit_lapl_egarch)[1]
    alpha11000_egarch_l[i] <- coef(fit_lapl_egarch)[2]
    beta11000_egarch_l[i] <- coef(fit_lapl_egarch)[3]
    gamma1000_egarch_l[i] <- coef(fit_lapl_egarch)[4]
    
    kounter=kounter+1
    i=i+1
    
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}

  if (i==pocetSim+1) {break}
  
}
print(j)

print(kounter)

print("omega")
print(summary(omega1000_egarch_l))
print(sd(omega1000_egarch_l))

print("alpha1")
print(summary(alpha11000_egarch_l))
print(sd(alpha11000_egarch_l))

print("beta1")
print(summary(beta11000_egarch_l))
print(sd(beta11000_egarch_l))

print("gamma")
print(summary(gamma1000_egarch_l))
print(sd(gamma1000_egarch_l))
###############################################################################################

###############  EGARCH(1,1) t 100 odhad parametru
###############################################################################################
spec1_t_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             distribution.model = "std",
                             fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                               "beta1" = 0.7,"gamma1"=0.28, "shape"=4))

show(spec1_t_egarch)

pocetSim=10000

j=100

omega100_egarch_t <- c()
alpha1100_egarch_t <- c()
beta1100_egarch_t <- c()
gamma100_egarch_t <-c()

kounter=0  
i=1

spec_t_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "std",
  fixed.pars = list("shape"=4))


repeat{
  x_obj_t_egarch <- ugarchpath(spec1_t_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_t_egarch <- x_obj_t_egarch@path$seriesSim
  
  try({
    
    fit_t_egarch <- ugarchfit(spec_t_egarch, data = x1_t_egarch)
    
    omega100_egarch_t[i] <- coef(fit_t_egarch)[1]
    alpha1100_egarch_t[i] <- coef(fit_t_egarch)[2]
    beta1100_egarch_t[i] <- coef(fit_t_egarch)[3]
    gamma100_egarch_t[i] <- coef(fit_t_egarch)[4]
    
    kounter=kounter+1
    i=i+1
    
    if (i%%1000==0){print(i)}
    
  } , silent = TRUE)
  if (i==pocetSim+1) {break}
  
}
print(j)

print(kounter)

print("omega")
print(summary(omega100_egarch_t))
print(sd(omega100_egarch_t))

print("alpha1")
print(summary(alpha1100_egarch_t))
print(sd(alpha1100_egarch_t))

print("beta1")
print(summary(beta1100_egarch_t))
print(sd(beta1100_egarch_t))

print("gamma")
print(summary(gamma100_egarch_t))
print(sd(gamma100_egarch_t))
###############################################################################################

###############  EGARCH(1,1) t 500 odhad parametru
###############################################################################################
spec1_t_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             distribution.model = "std",
                             fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                               "beta1" = 0.7,"gamma1"=0.28, "shape"=4))

show(spec1_t_egarch)

pocetSim=10000

j=500

omega500_egarch_t <- c()
alpha1500_egarch_t <- c()
beta1500_egarch_t <- c()
gamma500_egarch_t <-c()

kounter=0  
i=1

spec_t_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "std",
  fixed.pars = list("shape"=4))


repeat{
  x_obj_t_egarch <- ugarchpath(spec1_t_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_t_egarch <- x_obj_t_egarch@path$seriesSim
  
  try({
    
    fit_t_egarch <- ugarchfit(spec_t_egarch, data = x1_t_egarch)
    
    omega500_egarch_t[i] <- coef(fit_t_egarch)[1]
    alpha1500_egarch_t[i] <- coef(fit_t_egarch)[2]
    beta1500_egarch_t[i] <- coef(fit_t_egarch)[3]
    gamma500_egarch_t[i] <- coef(fit_t_egarch)[4]
    
    kounter=kounter+1
    i=i+1
    
    if (i%%1000==0){print(i)}
    
  } , silent = TRUE)
  
  if (i==pocetSim+1) {break}
  
}
print(j)

print(kounter)

print("omega")
print(summary(omega500_egarch_t))
print(sd(omega500_egarch_t))

print("alpha1")
print(summary(alpha1500_egarch_t))
print(sd(alpha1500_egarch_t))

print("beta1")
print(summary(beta1500_egarch_t))
print(sd(beta1500_egarch_t))

print("gamma")
print(summary(gamma500_egarch_t))
print(sd(gamma500_egarch_t))
###############################################################################################

###############  EGARCH(1,1) t 1000 odhad parametru
###############################################################################################
spec1_t_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             distribution.model = "std",
                             fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                               "beta1" = 0.7,"gamma1"=0.28, "shape"=4))

show(spec1_t_egarch)

pocetSim=10000

j=1000

omega1000_egarch_t <- c()
alpha11000_egarch_t <- c()
beta11000_egarch_t <- c()
gamma1000_egarch_t <-c()

kounter=0  
i=1

spec_t_egarch <- ugarchspec(variance.model = list(model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "std",
  fixed.pars = list("shape"=4))


repeat{
  x_obj_t_egarch <- ugarchpath(spec1_t_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_t_egarch <- x_obj_t_egarch@path$seriesSim
  
  try({
    
    fit_t_egarch <- ugarchfit(spec_t_egarch, data = x1_t_egarch)
    
    omega1000_egarch_t[i] <- coef(fit_t_egarch)[1]
    alpha11000_egarch_t[i] <- coef(fit_t_egarch)[2]
    beta11000_egarch_t[i] <- coef(fit_t_egarch)[3]
    gamma1000_egarch_t[i] <- coef(fit_t_egarch)[4]
    
    kounter=kounter+1
    i=i+1
    
    if (i%%1000==0){print(i)}
    
  } , silent = TRUE)

  if (i==pocetSim+1) {break}
  
}
print(j)

print(kounter)

print("omega")
print(summary(omega1000_egarch_t))
print(sd(omega1000_egarch_t))

print("alpha1")
print(summary(alpha11000_egarch_t))
print(sd(alpha11000_egarch_t))

print("beta1")
print(summary(beta11000_egarch_t))
print(sd(beta11000_egarch_t))

print("gamma")
print(summary(gamma1000_egarch_t))
print(sd(gamma1000_egarch_t))
###############################################################################################

###############  EGARCH(1,1) Normal 100 BT
###############################################################################################
spec1norm_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                               fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                 "beta1" = 0.7,"gamma1"=0.28))

show(spec1norm_egarch)

pocetSim=10000

j=100

zamitSBT_egarch = 0
zamitNSBT_egarch = 0
zamitPSBT_egarch = 0
zamitJoint_egarch = 0

kounter=0
i=1

spec_norm_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE))

repeat {
  
  x_obj_norm_egarch <- ugarchpath(spec1norm_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm_egarch <- x_obj_norm_egarch@path$seriesSim
  
  
  try({
    
    fit_norm_egarch <- ugarchfit(spec_norm_egarch, data = x1norm_egarch)
    
    #####Sign bias tests
    
    if(is(fit_norm_egarch, "uGARCHfilter")) z = fit_norm_egarch@filter$z else z = z = fit_norm_egarch@fit$z
    res = as.numeric(residuals(fit_norm_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_norm_egarch)

    SBTp_egarch <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch <- signb$prob[4]
    
    if (SBTp_egarch<0.05) {zamitSBT_egarch=zamitSBT_egarch+1}
    if (NSBTp_egarch<0.05) {zamitNSBT_egarch=zamitNSBT_egarch+1}
    if (PSBTp_egarch<0.05) {zamitPSBT_egarch=zamitPSBT_egarch+1}
    if (JointTp_egarch<0.05) {zamitJoint_egarch=zamitJoint_egarch+1}
    
    kounter=kounter+1

    i=i+1

    if (i%%1000==0){print(i)}
    
  } , silent = TRUE)
  
  if (i==pocetSim+1) {break}

}

ratioSBT_egarch=zamitSBT_egarch/pocetSim*100
ratioNSBT_egarch=zamitNSBT_egarch/pocetSim*100
ratioPSBT_egarch=zamitPSBT_egarch/pocetSim*100
ratioJoint_egarch=zamitJoint_egarch/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch <- c(ratioSBT_egarch,ratioNSBT_egarch,ratioPSBT_egarch,ratioJoint_egarch)

ratiomatice_egarch <- rbind(testy,ratio_egarch)

print(j)
print(ratiomatice_egarch)

print(kounter)

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

###############  EGARCH(1,1) Normal 500 BT
###############################################################################################
spec1norm_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                               fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                 "beta1" = 0.7,"gamma1"=0.28))

show(spec1norm_egarch)

pocetSim=10000

j=500

zamitSBT_egarch = 0
zamitNSBT_egarch = 0
zamitPSBT_egarch = 0
zamitJoint_egarch = 0

kounter=0
i=1

spec_norm_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE))

repeat {
  
  x_obj_norm_egarch <- ugarchpath(spec1norm_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm_egarch <- x_obj_norm_egarch@path$seriesSim
  
  
  try({
    
    fit_norm_egarch <- ugarchfit(spec_norm_egarch, data = x1norm_egarch)

    #####Sign bias tests
    
    if(is(fit_norm_egarch, "uGARCHfilter")) z = fit_norm_egarch@filter$z else z = z = fit_norm_egarch@fit$z
    res = as.numeric(residuals(fit_norm_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_norm_egarch)

    SBTp_egarch <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch <- signb$prob[4]
    
    if (SBTp_egarch<0.05) {zamitSBT_egarch=zamitSBT_egarch+1}
    if (NSBTp_egarch<0.05) {zamitNSBT_egarch=zamitNSBT_egarch+1}
    if (PSBTp_egarch<0.05) {zamitPSBT_egarch=zamitPSBT_egarch+1}
    if (JointTp_egarch<0.05) {zamitJoint_egarch=zamitJoint_egarch+1}
    
    kounter=kounter+1

    i=i+1

  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
}

ratioSBT_egarch=zamitSBT_egarch/pocetSim*100
ratioNSBT_egarch=zamitNSBT_egarch/pocetSim*100
ratioPSBT_egarch=zamitPSBT_egarch/pocetSim*100
ratioJoint_egarch=zamitJoint_egarch/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch <- c(ratioSBT_egarch,ratioNSBT_egarch,ratioPSBT_egarch,ratioJoint_egarch)

ratiomatice_egarch <- rbind(testy,ratio_egarch)

print(j)
print(ratiomatice_egarch)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) Normal 1000 BT
###############################################################################################
spec1norm_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                               fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                 "beta1" = 0.7,"gamma1"=0.28))

show(spec1norm_egarch)

pocetSim=10000

j=1000

zamitSBT_egarch = 0
zamitNSBT_egarch = 0
zamitPSBT_egarch = 0
zamitJoint_egarch = 0

kounter=0
i=1

spec_norm_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE))

repeat {
  
  x_obj_norm_egarch <- ugarchpath(spec1norm_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1norm_egarch <- x_obj_norm_egarch@path$seriesSim
  
  
  try({
    
    fit_norm_egarch <- ugarchfit(spec_norm_egarch, data = x1norm_egarch)

    #####Sign bias tests
    
    if(is(fit_norm_egarch, "uGARCHfilter")) z = fit_norm_egarch@filter$z else z = z = fit_norm_egarch@fit$z
    res = as.numeric(residuals(fit_norm_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_norm_egarch)

    SBTp_egarch <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch <- signb$prob[4]
    
    if (SBTp_egarch<0.05) {zamitSBT_egarch=zamitSBT_egarch+1}
    if (NSBTp_egarch<0.05) {zamitNSBT_egarch=zamitNSBT_egarch+1}
    if (PSBTp_egarch<0.05) {zamitPSBT_egarch=zamitPSBT_egarch+1}
    if (JointTp_egarch<0.05) {zamitJoint_egarch=zamitJoint_egarch+1}
    
    kounter=kounter+1

    i=i+1

  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}

}

ratioSBT_egarch=zamitSBT_egarch/pocetSim*100
ratioNSBT_egarch=zamitNSBT_egarch/pocetSim*100
ratioPSBT_egarch=zamitPSBT_egarch/pocetSim*100
ratioJoint_egarch=zamitJoint_egarch/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch <- c(ratioSBT_egarch,ratioNSBT_egarch,ratioPSBT_egarch,ratioJoint_egarch)

ratiomatice_egarch <- rbind(testy,ratio_egarch)

print(j)
print(ratiomatice_egarch)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) Laplace 100 BT
###############################################################################################
spec1_lapl_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                    garchOrder=c(1,1)),
                                mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                                distribution.model = "ged",
                                fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                  "beta1" = 0.7,"gamma1"=0.28, "shape"=1))

show(spec1_lapl_egarch)

pocetSim=10000

j=100

zamitSBT_egarch_lapl = 0
zamitNSBT_egarch_lapl = 0
zamitPSBT_egarch_lapl = 0
zamitJoint_egarch_lapl = 0

kounter=0  
i=1

spec_lapl_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "ged",
  fixed.pars = list("shape"=1))

repeat{
  
  x_obj_lapl_egarch <- ugarchpath(spec1_lapl_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl_egarch <- x_obj_lapl_egarch@path$seriesSim
  
  try({
    
    fit_lapl_egarch <- ugarchfit(spec_lapl_egarch, data = x1_lapl_egarch)

    #####Sign bias tests
    
    if(is(fit_lapl_egarch, "uGARCHfilter")) z = fit_lapl_egarch@filter$z else z = z = fit_lapl_egarch@fit$z
    res = as.numeric(residuals(fit_lapl_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_lapl_egarch)

    SBTp_egarch_lapl <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch_lapl <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch_lapl <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch_lapl <- signb$prob[4]
    
    if (SBTp_egarch_lapl<0.05) {zamitSBT_egarch_lapl=zamitSBT_egarch_lapl+1}
    if (NSBTp_egarch_lapl<0.05) {zamitNSBT_egarch_lapl=zamitNSBT_egarch_lapl+1}
    if (PSBTp_egarch_lapl<0.05) {zamitPSBT_egarch_lapl=zamitPSBT_egarch_lapl+1}
    if (JointTp_egarch_lapl<0.05) {zamitJoint_egarch_lapl=zamitJoint_egarch_lapl+1}
    
    kounter=kounter+1
    i=i+1

  } , silent = TRUE)

  if (i%%1000==0){print(i)}
  
  if (i==pocetSim+1) {break}
  
}

ratioSBT_egarch_lapl=zamitSBT_egarch_lapl/pocetSim*100
ratioNSBT_egarch_lapl=zamitNSBT_egarch_lapl/pocetSim*100
ratioPSBT_egarch_lapl=zamitPSBT_egarch_lapl/pocetSim*100
ratioJoint_egarch_lapl=zamitJoint_egarch_lapl/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch_lapl <- c(ratioSBT_egarch_lapl, ratioNSBT_egarch_lapl, ratioPSBT_egarch_lapl, ratioJoint_egarch_lapl)

ratiomatice_egarch_lapl <- rbind(testy,ratio_egarch_lapl)

print(j)
print(ratiomatice_egarch_lapl)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) Laplace 500 BT
###############################################################################################
spec1_lapl_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                    garchOrder=c(1,1)),
                                mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                                distribution.model = "ged",
                                fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                  "beta1" = 0.7,"gamma1"=0.28, "shape"=1))

show(spec1_lapl_egarch)

pocetSim=10000

j=500

zamitSBT_egarch_lapl = 0
zamitNSBT_egarch_lapl = 0
zamitPSBT_egarch_lapl = 0
zamitJoint_egarch_lapl = 0

kounter=0  
i=1

spec_lapl_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "ged",
  fixed.pars = list("shape"=1))


repeat{
  x_obj_lapl_egarch <- ugarchpath(spec1_lapl_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl_egarch <- x_obj_lapl_egarch@path$seriesSim
  
  try({
    
    fit_lapl_egarch <- ugarchfit(spec_lapl_egarch, data = x1_lapl_egarch)

    #####Sign bias tests
    
    if(is(fit_lapl_egarch, "uGARCHfilter")) z = fit_lapl_egarch@filter$z else z = z = fit_lapl_egarch@fit$z
    res = as.numeric(residuals(fit_lapl_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_lapl_egarch)

    SBTp_egarch_lapl <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch_lapl <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch_lapl <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch_lapl <- signb$prob[4]
    
    if (SBTp_egarch_lapl<0.05) {zamitSBT_egarch_lapl=zamitSBT_egarch_lapl+1}
    if (NSBTp_egarch_lapl<0.05) {zamitNSBT_egarch_lapl=zamitNSBT_egarch_lapl+1}
    if (PSBTp_egarch_lapl<0.05) {zamitPSBT_egarch_lapl=zamitPSBT_egarch_lapl+1}
    if (JointTp_egarch_lapl<0.05) {zamitJoint_egarch_lapl=zamitJoint_egarch_lapl+1}
    
    kounter=kounter+1
    i=i+1
    
  } , silent = TRUE)
  
  if (i%%1000==0){print(i)}

    if (i==pocetSim+1) {break}
  
}

ratioSBT_egarch_lapl=zamitSBT_egarch_lapl/pocetSim*100
ratioNSBT_egarch_lapl=zamitNSBT_egarch_lapl/pocetSim*100
ratioPSBT_egarch_lapl=zamitPSBT_egarch_lapl/pocetSim*100
ratioJoint_egarch_lapl=zamitJoint_egarch_lapl/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch_lapl <- c(ratioSBT_egarch_lapl, ratioNSBT_egarch_lapl, ratioPSBT_egarch_lapl, ratioJoint_egarch_lapl)

ratiomatice_egarch_lapl <- rbind(testy,ratio_egarch_lapl)

print(j)
print(ratiomatice_egarch_lapl)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) Laplace 1000 BT
###############################################################################################
spec1_lapl_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                    garchOrder=c(1,1)),
                                mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                                distribution.model = "ged",
                                fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                                  "beta1" = 0.7,"gamma1"=0.28, "shape"=1))

show(spec1_lapl_egarch)

pocetSim=10000

j=1000

zamitSBT_egarch_lapl = 0
zamitNSBT_egarch_lapl = 0
zamitPSBT_egarch_lapl = 0
zamitJoint_egarch_lapl = 0

kounter=0  
i=1

spec_lapl_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "ged",
  fixed.pars = list("shape"=1))


repeat{
  x_obj_lapl_egarch <- ugarchpath(spec1_lapl_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_lapl_egarch <- x_obj_lapl_egarch@path$seriesSim
  
  try({
    
    fit_lapl_egarch <- ugarchfit(spec_lapl_egarch, data = x1_lapl_egarch)
    
    #####Sign bias tests
    
    if(is(fit_lapl_egarch, "uGARCHfilter")) z = fit_lapl_egarch@filter$z else z = z = fit_lapl_egarch@fit$z
    res = as.numeric(residuals(fit_lapl_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_lapl_egarch)

    SBTp_egarch_lapl <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch_lapl <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch_lapl <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch_lapl <- signb$prob[4]
    
    
    if (SBTp_egarch_lapl<0.05) {zamitSBT_egarch_lapl=zamitSBT_egarch_lapl+1}
    if (NSBTp_egarch_lapl<0.05) {zamitNSBT_egarch_lapl=zamitNSBT_egarch_lapl+1}
    if (PSBTp_egarch_lapl<0.05) {zamitPSBT_egarch_lapl=zamitPSBT_egarch_lapl+1}
    if (JointTp_egarch_lapl<0.05) {zamitJoint_egarch_lapl=zamitJoint_egarch_lapl+1}
    
    kounter=kounter+1
    i=i+1
    
    
  } , silent = TRUE)
  
  if (i%%2000==0){print(i)}

  if (i==pocetSim+1) {break}
  
}

ratioSBT_egarch_lapl=zamitSBT_egarch_lapl/pocetSim*100
ratioNSBT_egarch_lapl=zamitNSBT_egarch_lapl/pocetSim*100
ratioPSBT_egarch_lapl=zamitPSBT_egarch_lapl/pocetSim*100
ratioJoint_egarch_lapl=zamitJoint_egarch_lapl/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch_lapl <- c(ratioSBT_egarch_lapl, ratioNSBT_egarch_lapl, ratioPSBT_egarch_lapl, ratioJoint_egarch_lapl)

ratiomatice_egarch_lapl <- rbind(testy,ratio_egarch_lapl)

print(j)
print(ratiomatice_egarch_lapl)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) t 100 BT
###############################################################################################
spec1_t_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             distribution.model = "std",
                             fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                               "beta1" = 0.7,"gamma1"=0.28, "shape"=4))

show(spec1_t_egarch)

pocetSim=10000

j=100

zamitSBT_egarch_t = 0
zamitNSBT_egarch_t = 0
zamitPSBT_egarch_t = 0
zamitJoint_egarch_t = 0

kounter=0  
i=1

spec_t_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "std",
  fixed.pars = list("shape"=4))


repeat{
  x_obj_t_egarch <- ugarchpath(spec1_t_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_t_egarch <- x_obj_t_egarch@path$seriesSim
  
  try({
    
    fit_t_egarch <- ugarchfit(spec_t_egarch, data = x1_t_egarch)

    #####Sign bias tests
    
    if(is(fit_t_egarch, "uGARCHfilter")) z = fit_t_egarch@filter$z else z = z = fit_t_egarch@fit$z
    res = as.numeric(residuals(fit_t_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_t_egarch)

    SBTp_egarch_t <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch_t <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch_t <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch_t <- signb$prob[4]
    
    if (SBTp_egarch_t<0.05) {zamitSBT_egarch_t=zamitSBT_egarch_t+1}
    if (NSBTp_egarch_t<0.05) {zamitNSBT_egarch_t=zamitNSBT_egarch_t+1}
    if (PSBTp_egarch_t<0.05) {zamitPSBT_egarch_t=zamitPSBT_egarch_t+1}
    if (JointTp_egarch_t<0.05) {zamitJoint_egarch_t=zamitJoint_egarch_t+1}
    
    kounter=kounter+1
    i=i+1
    
    if (i%%200==0){print(i)}
    
  } , silent = TRUE)

  if (i==pocetSim+1) {break}
  
}

ratioSBT_egarch_t=zamitSBT_egarch_t/pocetSim*100
ratioNSBT_egarch_t=zamitNSBT_egarch_t/pocetSim*100
ratioPSBT_egarch_t=zamitPSBT_egarch_t/pocetSim*100
ratioJoint_egarch_t=zamitJoint_egarch_t/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch_t <- c(ratioSBT_egarch_t, ratioNSBT_egarch_t, ratioPSBT_egarch_t, ratioJoint_egarch_t)

ratiomatice_egarch_t <- rbind(testy,ratio_egarch_t)

print(j)
print(ratiomatice_egarch_t)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) t 500 BT
###############################################################################################
spec1_t_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             distribution.model = "std",
                             fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                               "beta1" = 0.7,"gamma1"=0.28, "shape"=4))

show(spec1_t_egarch)

pocetSim=10000

j=500

zamitSBT_egarch_t = 0
zamitNSBT_egarch_t = 0
zamitPSBT_egarch_t = 0
zamitJoint_egarch_t = 0

kounter=0  
i=1

spec_t_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "std",
  fixed.pars = list("shape"=4))


repeat{
  x_obj_t_egarch <- ugarchpath(spec1_t_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_t_egarch <- x_obj_t_egarch@path$seriesSim
  
  try({
    
    fit_t_egarch <- ugarchfit(spec_t_egarch, data = x1_t_egarch)

    #####Sign bias tests
    
    if(is(fit_t_egarch, "uGARCHfilter")) z = fit_t_egarch@filter$z else z = z = fit_t_egarch@fit$z
    res = as.numeric(residuals(fit_t_egarch))
    z2 = z^2
    n = length(z)

    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_t_egarch)

    SBTp_egarch_t <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch_t <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch_t <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch_t <- signb$prob[4]
    
    if (SBTp_egarch_t<0.05) {zamitSBT_egarch_t=zamitSBT_egarch_t+1}
    if (NSBTp_egarch_t<0.05) {zamitNSBT_egarch_t=zamitNSBT_egarch_t+1}
    if (PSBTp_egarch_t<0.05) {zamitPSBT_egarch_t=zamitPSBT_egarch_t+1}
    if (JointTp_egarch_t<0.05) {zamitJoint_egarch_t=zamitJoint_egarch_t+1}
    
    kounter=kounter+1
    i=i+1
    
    if (i%%200==0){print(i)}
    
  } , silent = TRUE)

  if (i==pocetSim+1) {break}
  
}

ratioSBT_egarch_t=zamitSBT_egarch_t/pocetSim*100
ratioNSBT_egarch_t=zamitNSBT_egarch_t/pocetSim*100
ratioPSBT_egarch_t=zamitPSBT_egarch_t/pocetSim*100
ratioJoint_egarch_t=zamitJoint_egarch_t/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch_t <- c(ratioSBT_egarch_t, ratioNSBT_egarch_t, ratioPSBT_egarch_t, ratioJoint_egarch_t)

ratiomatice_egarch_t <- rbind(testy,ratio_egarch_t)

print(j)
print(ratiomatice_egarch_t)

print(kounter)
###############################################################################################

###############  EGARCH(1,1) t 1000 BT
###############################################################################################
spec1_t_egarch <- ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
                             distribution.model = "std",
                             fixed.pars = list("omega" = -0.15, "alpha1" = -0.23,
                                               "beta1" = 0.7,"gamma1"=0.28, "shape"=4))

show(spec1_t_egarch)

pocetSim=10000

j=1000

zamitSBT_egarch_t = 0
zamitNSBT_egarch_t = 0
zamitPSBT_egarch_t = 0
zamitJoint_egarch_t = 0

kounter=0  
i=1

spec_t_egarch <- ugarchspec(variance.model = list(#model="eGARCH",
  garchOrder=c(1,1)),
  mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
  distribution.model = "std",
  fixed.pars = list("shape"=4))


repeat{
  x_obj_t_egarch <- ugarchpath(spec1_t_egarch, n.sim = j, n.start = 1000, m.sim=1)
  
  x1_t_egarch <- x_obj_t_egarch@path$seriesSim
  
  try({
    
    fit_t_egarch <- ugarchfit(spec_t_egarch, data = x1_t_egarch)

    #####Sign bias tests
    
    if(is(fit_t_egarch, "uGARCHfilter")) z = fit_t_egarch@filter$z else z = z = fit_t_egarch@fit$z
    res = as.numeric(residuals(fit_t_egarch))
    z2 = z^2
    n = length(z)
    
    zminus = as.integer(res<0)
    zplus = 1-zminus
    czminus = zminus*res
    czplus = zplus*res
    cz = cbind(rep(1, n), zminus, czminus, czplus)
    cz = cz[1:(n-1),]
    z2 = matrix(z2[2:n],ncol = 1)
    cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
    
    fitSB = lm(y~1+zminus, data = cm)

    fitNSB = lm(y~1+czminus, data = cm)

    fitPSB = lm(y~1+czplus, data = cm)

    signb <- signbias(fit_t_egarch)

    SBTp_egarch_t <- summary(fitSB)$coefficients[2,4]
    NSBTp_egarch_t <- summary(fitNSB)$coefficients[2,4]
    PSBTp_egarch_t <- summary(fitPSB)$coefficients[2,4]
    JointTp_egarch_t <- signb$prob[4]
    
    if (SBTp_egarch_t<0.05) {zamitSBT_egarch_t=zamitSBT_egarch_t+1}
    if (NSBTp_egarch_t<0.05) {zamitNSBT_egarch_t=zamitNSBT_egarch_t+1}
    if (PSBTp_egarch_t<0.05) {zamitPSBT_egarch_t=zamitPSBT_egarch_t+1}
    if (JointTp_egarch_t<0.05) {zamitJoint_egarch_t=zamitJoint_egarch_t+1}
    
    kounter=kounter+1
    i=i+1
    
    if (i%%1000==0){print(i)}
    
  } , silent = TRUE)

  if (i==pocetSim+1) {break}
  
}

ratioSBT_egarch_t=zamitSBT_egarch_t/pocetSim*100
ratioNSBT_egarch_t=zamitNSBT_egarch_t/pocetSim*100
ratioPSBT_egarch_t=zamitPSBT_egarch_t/pocetSim*100
ratioJoint_egarch_t=zamitJoint_egarch_t/pocetSim*100

testy <- c("SBT","NSBT","PSBT","Joint")
ratio_egarch_t <- c(ratioSBT_egarch_t, ratioNSBT_egarch_t, ratioPSBT_egarch_t, ratioJoint_egarch_t)

ratiomatice_egarch_t <- rbind(testy,ratio_egarch_t)

print(j)
print(ratiomatice_egarch_t)

print(kounter)
###############################################################################################

################################################################################################
#####                                Realna data

getSymbols("MSFT", from ="2000-01-03", to = "2019-01-03")
colnames(MSFT)
start(MSFT)
end(MSFT)

MSFT = MSFT[, "MSFT.Adjusted", drop=F]

##### MSFT obrazek

x0=index(MSFT)
y0=as.numeric(MSFT)
df <- data.frame(MSFT)

ggplot(df, aes(x=x0,y=y0)) + geom_line() + theme_minimal() + ylab("") + xlab("") +
  ggtitle("MSFT") + scale_x_date(date_breaks = "2 year", date_labels = "%m-%Y") +
  scale_y_continuous(breaks = round(seq(20, 100, by = 20),1))#+

#####

# calculate log-returns for GARCH analysis
MSFT.ret = CalculateReturns(MSFT, method="log")

# remove first NA observation
MSFT.ret = MSFT.ret[-1,]

colnames(MSFT.ret) ="MSFT"

##### MSFT returns obrazek

x1=index(MSFT.ret)
y1=as.numeric(MSFT.ret)
df <- data.frame(MSFT.ret)

ggplot(df, aes(x=x1,y=y1)) + geom_line() + theme_minimal() + ylab("") + xlab("") +
  ggtitle("MSFT returns") + scale_x_date(date_breaks = "2 year", date_labels = "%m-%Y") +
  scale_y_continuous(breaks = round(seq(-0.15, 0.15, by = 0.05),3))#+

#####
##### Popisne statistiky, ADF test, Jarque-Bera test
length(MSFT.ret)

summary(MSFT.ret)
skewness(MSFT.ret)
kurtosis(MSFT.ret)
sd(MSFT.ret)

adf.test(MSFT.ret)

jarque.bera.test(MSFT.ret)

#####
##### Histogram MSFT returns
df <- data.frame(MSFT.ret)
ggplot(df, aes(x=MSFT.ret)) + 
  geom_histogram(binwidth=0.002, color="black", fill="white") + theme_minimal() + xlab("") +
  stat_function(fun = function(x) 
    dnorm(x, mean = mean(MSFT.ret), sd = sd(MSFT.ret)) * 0.002 * length(MSFT.ret), 
    color = "darkred", size = 1)

#####
##### Diagnosticke testy
garch11.spec = ugarchspec(variance.model = list(garchOrder=c(1,1)), 
                          mean.model = list(armaOrder=c(0,0)),
                          fixed.pars = list("mu"=0))

MSFT.garch11.fit = ugarchfit(spec=garch11.spec, data=MSFT.ret)

res=MSFT.garch11.fit@fit$residuals

z2=res^2
n = length(z2)
zminus = as.integer(res<0)
zplus = 1-zminus
czminus = zminus*res
czplus = zplus*res
cz = cbind(rep(1, n), zminus, czminus, czplus)
cz = cz[1:(n-1),]
z2 = matrix(z2[2:n],ncol = 1)
cm = data.frame(y = z2, const = cz[,1], zminus = cz[,2], czminus = cz[,3], czplus = cz[,4])
fitA = lm(y~1+zminus+czminus+czplus, data = cm)
resA = residuals(fitA)
sgmatrix = matrix(ncol = 3, nrow = 4)
rownames(sgmatrix) = c("Sign Bias","Negative Sign Bias","Positive Sign Bias","Joint Effect")
colnames(sgmatrix) = c("t-value","prob","sig")
sgmatrix[1:3,1:2] = abs(summary(fitA)$coef[2:4,3:4])
jeffect = linearHypothesis(fitA, c("zminus = 0", "czminus = 0","czplus  =  0"), test  =  "Chisq")
sgmatrix[4,1] = jeffect[5]$Chisq[2]
sgmatrix[4,2] = jeffect[6][2,]
sgmatrix = as.data.frame(sgmatrix)

print(sgmatrix)

fitSB = lm(y~1+zminus, data = cm)
summary(fitSB)

fitNSB = lm(y~1+czminus, data = cm)
summary(fitNSB)

fitPSB = lm(y~1+czplus, data = cm)
summary(fitPSB)

#####
##### GARCH(1,1) Normal distr
garch11.spec = ugarchspec(variance.model = list(garchOrder=c(1,1)), 
                          mean.model = list(armaOrder=c(0,0)),
                          fixed.pars = list("mu"=0))

MSFT.garch11.fit = ugarchfit(spec=garch11.spec, data=MSFT.ret)
MSFT.garch11.fit

##### EGARCH(1,1) Normal distr
egarch11.spec = ugarchspec(variance.model=list(model="eGARCH",
                                               garchOrder=c(1,1)),
                           mean.model=list(armaOrder=c(0,0)),
                           fixed.pars = list("mu"=0))

MSFT.egarch11.fit = ugarchfit(egarch11.spec, MSFT.ret)
MSFT.egarch11.fit

##### NIC GARCH EGARCH Normal distr
nic.garch11 = newsimpact(MSFT.garch11.fit)
nic.egarch11 = newsimpact(MSFT.egarch11.fit)

x  <- nic.egarch11$zx
y1 <- nic.egarch11$zy
y2 <- nic.garch11$zy
df <- data.frame(x,y1,y2)
xlabel<-expression(epsilon[t])

ggplot(df, aes(x)) +                    
  geom_line(aes(y=y1, linetype = "EGARCH(1,1)"), colour="black") +  
  geom_line(aes(y=y2, linetype = "GARCH(1,1)"), colour="black") + 
  geom_vline(xintercept = 0) + 
  ylab("NIC") + xlab(xlabel) + 
  theme_classic() +  theme( legend.title = element_blank(),
                            axis.text.y = element_blank()
  ) +
  scale_colour_manual(name="", values=c("EGARCH(1,1)"='solid', "GARCH(1,1)"='dashed')) 

##### GARCH(1,1)  with Student-t errors
garch11.t.spec = ugarchspec(variance.model = list(garchOrder=c(1,1)), 
                            mean.model = list(armaOrder=c(0,0)),
                            distribution.model = "std",
                            fixed.pars = list("mu"=0,"shape"=4 ))

MSFT.garch11.t.fit = ugarchfit(spec=garch11.t.spec, data=MSFT.ret)                             
MSFT.garch11.t.fit

##### EGARCH(1,1)  with Student-t errors
egarch11.spec.t = ugarchspec(variance.model=list(model="eGARCH",
                                                 garchOrder=c(1,1)),
                             mean.model=list(armaOrder=c(0,0)),
                             distribution.model = "std",
                             fixed.pars = list("mu"=0,"shape"=4))

MSFT.egarch11.t.fit = ugarchfit(egarch11.spec.t, MSFT.ret)
MSFT.egarch11.t.fit

##### NIC GARCH EGARCH t distr

nic.garch11.t.fit = newsimpact(MSFT.garch11.t.fit)
nic.egarch11.t.fit = newsimpact(MSFT.egarch11.t.fit)

x.t  <- nic.egarch11.t.fit$zx
y1.t <- nic.egarch11.t.fit$zy
y2.t <- nic.garch11.t.fit$zy
df.t <- data.frame(x.t,y1.t,y2.t)
model<- c("GARCH(1,1)","EGARCH(1,1)")
xlabel<-expression(epsilon[t])

ggplot(df.t, aes(x=x.t)) +                    
  geom_line(aes(y=y1.t, linetype = "EGARCH(1,1)"), colour="black") +  
  geom_line(aes(y=y2.t, linetype = "GARCH(1,1)"), colour="black") + 
  geom_vline(xintercept = 0) + 
  ylab("NIC") + xlab(xlabel) + 
  theme_classic() +  theme( legend.title = element_blank(),
                            axis.text.y = element_blank()) +
  scale_colour_manual(name="", values=c("EGARCH(1,1)"='solid', "GARCH(1,1)"='dashed'))

##### GARCH(1,1)  with GED errors
garch11.ged.spec = ugarchspec(variance.model = list(garchOrder=c(1,1)), 
                              mean.model = list(armaOrder=c(0,0)),
                              distribution.model = "ged",
                              fixed.pars = list("shape" = 1,"mu"=0))

MSFT.garch11.ged.fit = ugarchfit(spec=garch11.ged.spec, data=MSFT.ret)                             
MSFT.garch11.ged.fit

##### EGARCH(1,1)  with GED errors
egarch11.spec.ged = ugarchspec(variance.model=list(model="eGARCH",
                                                   garchOrder=c(1,1)),
                               mean.model=list(armaOrder=c(0,0)),
                               distribution.model = "ged",
                               fixed.pars = list("shape" = 1,"mu"=0))

MSFT.egarch11.ged.fit = ugarchfit(egarch11.spec.ged, MSFT.ret)
MSFT.egarch11.ged.fit

##### NIC GARCH EGARCH Laplace distr
nic.garch11.ged.fit = newsimpact(MSFT.garch11.ged.fit)
nic.egarch11.ged.fit = newsimpact(MSFT.egarch11.ged.fit)

x.ged  <- nic.egarch11.ged.fit$zx
y1.ged <- nic.egarch11.ged.fit$zy
y2.ged <- nic.garch11.ged.fit$zy
df.ged <- data.frame(x.ged,y1.ged,y2.ged)
model<- c("GARCH(1,1)","EGARCH(1,1)")
xlabel<-expression(epsilon[t])

ggplot(df.ged, aes(x=x.ged)) +                    
  geom_line(aes(y=y1.ged, linetype = "EGARCH(1,1)"), colour="black") +  
  geom_line(aes(y=y2.ged, linetype = "GARCH(1,1)"), colour="black") + 
  geom_vline(xintercept = 0) + 
  ylab("NIC") + xlab(xlabel) + 
  theme_classic() +  theme( legend.title = element_blank(),
                            axis.text.y = element_blank()) +
  scale_colour_manual(name="", values=c("EGARCH(1,1)"='solid', "GARCH(1,1)"='dashed'))