0. LOAD PACKAGE & DATA
library(rugarch); library(fGarch); library(readxl)
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
##
## sigma
## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")
Data.xlsx <- read_excel("C:/Users/84896/Desktop/Tutorial_VaR/Data.xlsx")
Data <- Data.xlsx[,2:ncol(Data.xlsx)]
Date <- as.Date(Data.xlsx$Date)
1. ARMA MODEL
arma.model <- function(x) {
fits <- lapply(x, function(series) {
model <- autoarfima(series, include.mean = TRUE, ar.max = 2, ma.max = 2, criterion = "BIC", method = "full")
coef <- model$fit@fit$coef
bic <- model$rank.matrix[1,7]
mu <- model$fit@fit$coef['mu']
ar1 <- if ("ar1" %in% names(coef)) coef["ar1"] else NA
ar2 <- if ("ar2" %in% names(coef)) coef["ar2"] else NA
ma1 <- if ("ma1" %in% names(coef)) coef["ma1"] else NA
ma2 <- if ("ma2" %in% names(coef)) coef["ma2"] else NA
arma.coef <- data.frame(Mean = mu, AR_1 = ar1, AR_2 = ar2, MA_1 = ma1, MA_2 = ma2, BIC = bic)
arma.res <- model$fit@fit$residuals
list(coef = arma.coef, residuals = arma.res)
})
names(fits) <- names(x)
return(fits)
}
arma.fit <- arma.model(Data)
arma.fit <- list(coef = setNames(lapply(arma.fit, `[[`, "coef"), names(arma.fit)), residuals = setNames(lapply(arma.fit, `[[`, "residuals"), names(arma.fit)))
data.frame(Variable = rownames(do.call(rbind, arma.fit$coef)), do.call(rbind, arma.fit$coef), row.names = NULL)
## Variable Mean AR_1 AR_2 MA_1 MA_2 BIC
## 1 BRO 0.0005656314 NA NA NA NA -5.395531
## 2 CMS 0.0001121981 NA NA NA NA -5.138206
## 3 F -0.0001426368 0 0.04641178 NA NA -4.416993
## 4 JBL 0.0002153790 NA NA NA NA -4.083370
## 5 COST 0.0004467956 NA NA NA NA -5.262070
2. GARCH MODEL
SelectGARCH <- function(x){
garch11.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(1,1), model="sGARCH"), distribution.model= "norm")
egarch11.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(1,1), model="eGARCH"), distribution.model= "norm")
tgarch11.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(1,1), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm")
garch12.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(1,2), model="sGARCH"), distribution.model= "norm")
egarch12.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(1,2), model="eGARCH"), distribution.model= "norm")
tgarch12.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(1,2), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm")
garch21.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(2,1), model="sGARCH"), distribution.model= "norm")
egarch21.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(2,1), model="eGARCH"), distribution.model= "norm")
tgarch21.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(2,1), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm")
garch22.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(2,2), model="sGARCH"), distribution.model= "norm")
egarch22.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(2,2), model="eGARCH"), distribution.model= "norm")
tgarch22.model <- ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE), variance.model= list(garchOrder=c(2,2), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm")
garch11.fit <- ugarchfit(spec = garch11.model, data = x)
egarch11.fit <- ugarchfit(spec = egarch11.model, data = x)
tgarch11.fit <- ugarchfit(spec = tgarch11.model, data = x)
garch12.fit <- ugarchfit(spec = garch12.model, data = x)
egarch12.fit <- ugarchfit(spec = egarch12.model, data = x)
tgarch12.fit <- ugarchfit(spec = tgarch12.model, data = x)
garch21.fit <- ugarchfit(spec = garch21.model, data = x)
egarch21.fit <- ugarchfit(spec = egarch21.model, data = x)
tgarch21.fit <- ugarchfit(spec = tgarch21.model, data = x)
garch22.fit <- ugarchfit(spec = garch22.model, data = x)
egarch22.fit <- ugarchfit(spec = egarch22.model, data = x)
tgarch22.fit <- ugarchfit(spec = tgarch22.model, data = x)
GARCH11 <- infocriteria(garch11.fit)[2]
EGARCH11 <- infocriteria(egarch11.fit)[2]
TGARCH11 <- infocriteria(tgarch11.fit)[2]
GARCH12 <- infocriteria(garch12.fit)[2]
EGARCH12 <- infocriteria(egarch12.fit)[2]
TGARCH12 <- infocriteria(tgarch12.fit)[2]
GARCH21 <- infocriteria(garch21.fit)[2]
EGARCH21 <- infocriteria(egarch21.fit)[2]
TGARCH21 <- infocriteria(tgarch21.fit)[2]
GARCH22 <- infocriteria(garch22.fit)[2]
EGARCH22 <- infocriteria(egarch22.fit)[2]
TGARCH22 <- infocriteria(tgarch22.fit)[2]
data.frame('GARCH(1,1)'=GARCH11, 'GARCH(1,2)'=GARCH12, 'GARCH(2,1)'=GARCH21, 'GARCH(2,2)'=GARCH22,
'EGARCH(1,1)'=EGARCH11, 'EGARCH(1,2)'=EGARCH12, 'EGARCH(2,1)'=EGARCH21, 'EGARCH(2,2)'=EGARCH22,
'TGARCH(1,1)'=TGARCH11, 'TGARCH(1,2)'=TGARCH12, 'TGARCH(2,1)'=TGARCH21, 'TGARCH(2,2)'=TGARCH22)
}
RunGARCHForAll <- function(data_list) {
results <- lapply(names(data_list), function(name) {
bic_values <- SelectGARCH(data_list[[name]])
bic_values$Variable <- name
return(bic_values)
})
final_results <- do.call(rbind, results)
final_results <- final_results[, c("Variable", setdiff(names(final_results), "Variable"))]
return(final_results)
}
df_IC <- as.data.frame(t(RunGARCHForAll(arma.fit$residuals)))
colnames(df_IC) <- df_IC[1, ]
model_IC <- df_IC[-1, ]
print(model_IC)
## BRO CMS F JBL COST
## GARCH.1.1. -5.654710 -5.684385 -4.768260 -4.437619 -5.575247
## GARCH.1.2. -5.654092 -5.683521 -4.771919 -4.437058 -5.573753
## GARCH.2.1. -5.653310 -5.683343 -4.766799 -4.436120 -5.572882
## GARCH.2.2. -5.652650 -5.682072 -4.770477 -4.435616 -5.572310
## EGARCH.1.1. -5.668553 -5.693706 -4.778308 -4.487293 -5.595543
## EGARCH.1.2. -5.670116 -5.692735 -4.782375 -4.485748 -5.593283
## EGARCH.2.1. -5.669208 -5.692396 -4.786770 -4.486374 -5.592522
## EGARCH.2.2. -5.669490 -5.691297 -4.785390 -4.484981 -5.592001
## TGARCH.1.1. -5.671972 -5.691121 -4.779967 -4.491415 -5.595696
## TGARCH.1.2. -5.673457 -5.689669 -4.783524 -4.489028 -5.591139
## TGARCH.2.1. -5.669342 -5.694353 -4.776568 -4.487864 -5.589518
## TGARCH.2.2. -5.671451 -5.696338 -4.780862 -4.486650 -5.590903
print(data.frame(Best_Model= apply(model_IC, 2, function(col) rownames(model_IC)[which.min(col)])))
## Best_Model
## BRO TGARCH.1.2.
## CMS TGARCH.2.2.
## F EGARCH.2.1.
## JBL TGARCH.1.1.
## COST TGARCH.1.1.
BRO.fit <- ugarchfit(data = arma.fit$residuals$BRO,
spec = ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE),
variance.model= list(garchOrder=c(1,2), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm"))
CMS.fit <- ugarchfit(data = arma.fit$residuals$CMS,
spec = ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE),
variance.model= list(garchOrder=c(2,2), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm"))
F.fit <- ugarchfit(data = arma.fit$residuals$'F',
spec = ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE),
variance.model= list(garchOrder=c(2,1), model="eGARCH"), distribution.model= "norm"))
JBL.fit <- ugarchfit(data = arma.fit$residuals$JBL,
spec = ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE),
variance.model= list(garchOrder=c(1,1), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm"))
COST.fit <- ugarchfit(data = arma.fit$residuals$COST,
spec = ugarchspec(mean.model= list(armaOrder= c(0,0), include.mean= FALSE),
variance.model= list(garchOrder=c(1,1), model="fGARCH", submodel= "TGARCH"), distribution.model= "norm"))
3. FORECAST DAILY VARIANCE
# COUNT TRADING DAYS
date_pred <- seq(as.Date("2024-01-01"), as.Date("2024-03-31"), by = "day")
trade_day <- length(date_pred[!(weekdays(date_pred) %in% c("Saturday", "Sunday"))])
date_pred <- date_pred[!(weekdays(date_pred) %in% c("Saturday", "Sunday"))]
merged_dates <- sort(union(Date, date_pred))
full_dates <- as.Date(merged_dates, origin = "1970-01-01")
# FITTED VARIANCE
variance_BRO <- BRO.fit@fit$var
variance_CMS <- CMS.fit@fit$var
variance_F <- F.fit@fit$var
variance_JBL <- JBL.fit@fit$var
variance_COST <- COST.fit@fit$var
# FORECASE DAILY VARIANCE
fspec_BRO <- getspec(BRO.fit) ### BRO ###
setfixed(fspec_BRO) <- as.list(coef(BRO.fit))
pred_BRO <- ugarchforecast(fspec_BRO, data = arma.fit$residuals$BRO, n.ahead = 1, n.roll = trade_day-1, out.sample = trade_day-1)
variance_pred_BRO <- pred_BRO@forecast$sigmaFor^2
fspec_CMS <- getspec(CMS.fit) ### CMS ###
setfixed(fspec_CMS) <- as.list(coef(CMS.fit))
pred_CMS <- ugarchforecast(fspec_CMS, data = arma.fit$residuals$CMS, n.ahead = 1, n.roll = trade_day-1, out.sample = trade_day-1)
variance_pred_CMS <- pred_CMS@forecast$sigmaFor^2
fspec_F <- getspec(F.fit) ### F ###
setfixed(fspec_F) <- as.list(coef(F.fit))
pred_F <- ugarchforecast(fspec_F, data = arma.fit$residuals$'F', n.ahead = 1, n.roll = trade_day-1, out.sample = trade_day-1)
variance_pred_F <- pred_F@forecast$sigmaFor^2
fspec_JBL <- getspec(JBL.fit) ### JBL ###
setfixed(fspec_JBL) <- as.list(coef(JBL.fit))
pred_JBL <- ugarchforecast(fspec_JBL, data = arma.fit$residuals$JBL, n.ahead = 1, n.roll = trade_day-1, out.sample = trade_day-1)
variance_pred_JBL <- pred_JBL@forecast$sigmaFor^2
fspec_COST <- getspec(COST.fit) ### COST ###
setfixed(fspec_COST) <- as.list(coef(COST.fit))
pred_COST <- ugarchforecast(fspec_COST, data = arma.fit$residuals$COST, n.ahead = 1, n.roll = trade_day-1, out.sample = trade_day-1)
variance_pred_COST <- pred_COST@forecast$sigmaFor^2
# DISPLAY RESULTS
plot(Date, variance_BRO, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "BRO")
lines(date_pred, variance_pred_BRO, col = "red", lty = 2) # Đường nét đứt cho dự báo
legend("topleft", legend = c("Fitted Variance", "Predicted Variance"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, variance_CMS, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "CMS")
lines(date_pred, variance_pred_CMS, col = "red", lty = 2) # Đường nét đứt cho dự báo
legend("topleft", legend = c("Fitted Variance", "Predicted Variance"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, variance_F, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "F")
lines(date_pred, variance_pred_F, col = "red", lty = 2) # Đường nét đứt cho dự báo
legend("topleft", legend = c("Fitted Variance", "Predicted Variance"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, variance_JBL, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "JBL")
lines(date_pred, variance_pred_JBL, col = "red", lty = 2) # Đường nét đứt cho dự báo
legend("topleft", legend = c("Fitted Variance", "Predicted Variance"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, variance_COST, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "COST")
lines(date_pred, variance_pred_COST, col = "red", lty = 2) # Đường nét đứt cho dự báo
legend("topleft", legend = c("Fitted Variance", "Predicted Variance"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
4. FORECAST DAILY VALUE AT RISK
VaR_BRO <- as.numeric(quantile(BRO.fit, probs = 0.01))
VaR_pred_BRO <- as.numeric(quantile(pred_BRO, probs = 0.01))
VaR_CMS <- as.numeric(quantile(CMS.fit, probs = 0.01))
VaR_pred_CMS <- as.numeric(quantile(pred_CMS, probs = 0.01))
VaR_F <- as.numeric(quantile(F.fit, probs = 0.01))
VaR_pred_F <- as.numeric(quantile(pred_F, probs = 0.01))
VaR_JBL <- as.numeric(quantile(JBL.fit, probs = 0.01))
VaR_pred_JBL <- as.numeric(quantile(pred_JBL, probs = 0.01))
VaR_COST <- as.numeric(quantile(COST.fit, probs = 0.01))
VaR_pred_COST <- as.numeric(quantile(pred_COST, probs = 0.01))
# DISPLAY RESULTS
plot(Date, VaR_BRO, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "BRO")
lines(date_pred, VaR_pred_BRO, col = "red", lty = 2)
legend("bottomleft", legend = c("Fitted VaR", "Predicted VaR"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, VaR_CMS, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "CMS")
lines(date_pred, VaR_pred_CMS, col = "red", lty = 2)
legend("bottomleft", legend = c("Fitted VaR", "Predicted VaR"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, VaR_F, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "F")
lines(date_pred, VaR_pred_F, col = "red", lty = 2)
legend("bottomleft", legend = c("Fitted VaR", "Predicted VaR"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, VaR_JBL, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "JBL")
lines(date_pred, VaR_pred_JBL, col = "red", lty = 2)
legend("bottomleft", legend = c("Fitted VaR", "Predicted VaR"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
plot(Date, VaR_COST, type = "l", col = "black", xlab = "Time", ylab = "Value at Risk", main = "COST")
lines(date_pred, VaR_pred_COST, col = "red", lty = 2)
legend("bottomleft", legend = c("Fitted VaR", "Predicted VaR"), col = c("black", "red"), lty = c(1, 2), lwd = 3, bty = "n")
5. GLOBAL MINIMUM VARIANCE PORTFOLIO