Libraries
library(dplyr)
library(psych)
library(readxl)
library(writexl)
library(kableExtra)
# Econometrics
library(tidyverse)
library(plm) # panel models
library(sandwich) #covariance matrixes
library(lmtest) # tests
library(xtable)# latex tables
library(stargazer) # latex regression tables
library(ggpubr) # correlation test
# for this HA
library(tsDyn)
library(tseries)
library(tidyverse)
library(PerformanceAnalytics)
library(rugarch)
library(forecast)
dat_1 <- read_excel(
'C:/Users/Popov/Documents/Studies/NES_studies/R/Time_Series/HA3/European_indexes.xls')
describe(dat_1)
## vars n mean sd median trimmed mad min max range
## ...1 1 2542 NaN NA NA NaN NA Inf -Inf -Inf
## DAX 2 2542 3320.03 1868.31 2308.20 3057.58 1043.90 1322.7 8136.16 6813.46
## CAC 3 2542 2903.82 1465.65 2118.60 2621.99 477.92 1441.0 6944.77 5503.77
## FTSE 4 2542 4140.38 1465.12 3662.55 4055.99 1579.64 2054.8 6950.60 4895.80
## skew kurtosis se
## ...1 NA NA NA
## DAX 0.91 -0.47 37.06
## CAC 1.34 0.55 29.07
## FTSE 0.48 -1.30 29.06
Data cleaning
dat_1 <- dat_1[, !names(dat_1) %in% "...1"] # null anyway
cac_ts <- ts(dat_1$CAC, frequency = 1)
variables creation
# returns
ret <- Return.calculate(cac_ts, method = "discrete")[2:2542]
# create train-test samples, sq roots
sqrt_abs_ret <- ts(sqrt(1e-3 + abs(ret)))
train_sqrt_abs <- ts(sqrt_abs_ret[1:1500])
test_sqrt_abs <- ts(sqrt_abs_ret[1501:2541])
# create train-test samples, log
log_abs_ret <- log(1e-3 + abs(ret))
train_log_abs <- ts(log_abs_ret[1:1500])
test_log_abs <- ts(log_abs_ret[1501:2541])
# create train-test samples, diff log
dlog_cac <- diff(log(1e-5 + cac_ts))
train_dlog <- ts(dlog_cac[1:1500])
test_dlog <- ts(dlog_cac[1501:2541])
Problem 1
# Model 1: sGARCH
g1 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = 'norm')
sGARCH_fit <- ugarchfit(spec = g1, data = train_dlog)
coef_sGARCH <- coef(sGARCH_fit)
robust_se_sGARCH <- sGARCH_fit@fit$robust.se.coef
t_stats_sGARCH <- coef_sGARCH / robust_se_sGARCH
# Model 2: gjrGARCH
g2 <- ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = 'norm')
gjrGARCH_fit <- ugarchfit(spec = g2, data = train_dlog)
coef_gjrGARCH <- coef(gjrGARCH_fit)
robust_se_gjrGARCH <- gjrGARCH_fit@fit$robust.se.coef
t_stats_gjrGARCH <- coef_gjrGARCH / robust_se_gjrGARCH
# Model 3: GARCH-in-mean
g3 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE, archm = TRUE))
garchM_fit <- ugarchfit(spec = g3, data = train_dlog)
coef_garchM <- coef(garchM_fit)
robust_se_garchM <- garchM_fit@fit$robust.se.coef
t_stats_garchM <- coef_garchM / robust_se_garchM
# DataFrame
results <- data.frame(
Model = rep(c("sGARCH", "gjrGARCH", "GARCH-in-mean"),
times = c(length(coef_sGARCH), length(coef_gjrGARCH), length(coef_garchM))),
Parameter = c(names(coef_sGARCH), names(coef_gjrGARCH), names(coef_garchM)),
Coefficient = c(coef_sGARCH, coef_gjrGARCH, coef_garchM),
Robust_SE = c(robust_se_sGARCH, robust_se_gjrGARCH, robust_se_garchM),
T_Statistic = c(t_stats_sGARCH, t_stats_gjrGARCH, t_stats_garchM)
)
# Results
print(results)
## Model Parameter Coefficient Robust_SE T_Statistic
## 1 sGARCH mu 2.886003e-04 2.588807e-04 1.1148001
## 2 sGARCH omega 7.137552e-06 3.519936e-07 20.2775049
## 3 sGARCH alpha1 7.202903e-02 3.277691e-03 21.9755432
## 4 sGARCH beta1 8.630308e-01 8.618321e-03 100.1390867
## 5 gjrGARCH mu 1.274115e-04 2.579391e-04 0.4939595
## 6 gjrGARCH omega 5.354480e-06 2.792636e-07 19.1735733
## 7 gjrGARCH alpha1 1.297214e-02 6.473438e-03 2.0039029
## 8 gjrGARCH beta1 8.953681e-01 7.017286e-03 127.5946430
## 9 gjrGARCH gamma1 8.599411e-02 2.215624e-02 3.8812588
## 10 GARCH-in-mean mu -9.185485e-03 3.671951e-03 -2.5015271
## 11 GARCH-in-mean archm 9.400415e-01 3.433525e-01 2.7378325
## 12 GARCH-in-mean omega 5.994246e-06 2.863273e-07 20.9349450
## 13 GARCH-in-mean alpha1 4.974871e-02 5.350657e-03 9.2976823
## 14 GARCH-in-mean beta1 8.932530e-01 9.552319e-03 93.5116365
Problem 2
# results dataframe
results_forecast <- data.frame(
Model = character(),
MSPE_m = numeric(),
MSPE_v = numeric(),
stringsAsFactors = FALSE
)
# sGARCH Forecast
g1f <- getspec(sGARCH_fit)
setfixed(g1f) <- as.list(coef(sGARCH_fit))
sGARCH_forecast <- ugarchforecast(g1f, n.ahead = 1, n.roll = 1040, data = test_dlog, out.sample = 1040)
sigma_sGARCH <- sGARCH_forecast@forecast$sigmaFor[1, 1:1041]
fitted_sGARCH <- sGARCH_forecast@forecast$seriesFor[1, 1:1041]
residuals_sGARCH <- test_dlog - fitted_sGARCH
MSPE_m_sGARCH <- sum(residuals_sGARCH^2) / length(residuals_sGARCH)
MSPE_v_sGARCH <- sum((residuals_sGARCH^2 - sigma_sGARCH^2)^2) / length(residuals_sGARCH)
# Add sGARCH results to the table
results_forecast <- rbind(results_forecast, data.frame(
Model = "sGARCH",
MSPE_m = MSPE_m_sGARCH,
MSPE_v = MSPE_v_sGARCH
))
# gjrGARCH Forecast
g2f <- getspec(gjrGARCH_fit)
setfixed(g2f) <- as.list(coef(gjrGARCH_fit))
gjrGARCH_forecast <- ugarchforecast(g2f, n.ahead = 1, n.roll = 1040, data = test_dlog, out.sample = 1040)
sigma_gjrGARCH <- gjrGARCH_forecast@forecast$sigmaFor[1, 1:1041]
fitted_gjrGARCH <- gjrGARCH_forecast@forecast$seriesFor[1, 1:1041]
residuals_gjrGARCH <- test_dlog - fitted_gjrGARCH
MSPE_m_gjrGARCH <- sum(residuals_gjrGARCH^2) / length(residuals_gjrGARCH)
MSPE_v_gjrGARCH <- sum((residuals_gjrGARCH^2 - sigma_gjrGARCH^2)^2) / length(residuals_gjrGARCH)
# Add gjrGARCH results to the table
results_forecast <- rbind(results_forecast, data.frame(
Model = "gjrGARCH",
MSPE_m = MSPE_m_gjrGARCH,
MSPE_v = MSPE_v_gjrGARCH
))
# GARCH-in-mean Forecast
g3f <- getspec(garchM_fit)
setfixed(g3f) <- as.list(coef(garchM_fit))
garchM_forecast <- ugarchforecast(g3f, n.ahead = 1, n.roll = 1040, data = test_dlog, out.sample = 1040)
sigma_garchM <- garchM_forecast@forecast$sigmaFor[1, 1:1041]
fitted_garchM <- garchM_forecast@forecast$seriesFor[1, 1:1041]
residuals_garchM <- test_dlog - fitted_garchM
MSPE_m_garchM <- sum(residuals_garchM^2) / length(residuals_garchM)
MSPE_v_garchM <- sum((residuals_garchM^2 - sigma_garchM^2)^2) / length(residuals_garchM)
# Add GARCH-in-mean results to the table
results_forecast <- rbind(results_forecast, data.frame(
Model = "GARCH-in-mean",
MSPE_m = MSPE_m_garchM,
MSPE_v = MSPE_v_garchM
))
# Print the results table
print(results_forecast)
## Model MSPE_m MSPE_v
## 1 sGARCH 0.0001903806 1.075183e-07
## 2 gjrGARCH 0.0001906156 1.066751e-07
## 3 GARCH-in-mean 0.0002016830 1.574548e-07
Problem 3
# sGARCH (GED) model
g4 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = 'ged')
garch11_ged <- ugarchfit(spec = g4, data = train_dlog)
# Extract results
coef_garch11_ged <- coef(garch11_ged)
robust_se_garch11_ged <- garch11_ged@fit$robust.se.coef
t_stats_garch11_ged <- coef_garch11_ged / robust_se_garch11_ged
# results table
results_ged <- data.frame(
Model = rep("sGARCH (GED)", length(coef_garch11_ged)),
Parameter = names(coef_garch11_ged),
Coefficient = coef_garch11_ged,
Robust_SE = robust_se_garch11_ged,
T_Statistic = t_stats_garch11_ged,
row.names = NULL # Remove row indices
)
print(results_ged, row.names = FALSE)
## Model Parameter Coefficient Robust_SE T_Statistic
## sGARCH (GED) mu 2.616525e-04 2.409505e-04 1.085918
## sGARCH (GED) omega 6.121878e-06 4.377929e-07 13.983504
## sGARCH (GED) alpha1 6.417396e-02 3.725314e-03 17.226456
## sGARCH (GED) beta1 8.802417e-01 7.592750e-03 115.931875
## sGARCH (GED) shape 1.592960e+00 9.481961e-02 16.799904
Problem 4
library(forecast)
# Fit AR models
Ar1 <- Arima(train_log_abs, order = c(1, 0, 0), include.mean = TRUE, method = "ML")
Ar2 <- Arima(train_log_abs, order = c(2, 0, 0), include.mean = TRUE, method = "ML")
Ar3 <- Arima(train_log_abs, order = c(3, 0, 0), include.mean = TRUE, method = "ML")
Ar4 <- Arima(train_log_abs, order = c(4, 0, 0), include.mean = TRUE, method = "ML")
Ar5 <- Arima(train_log_abs, order = c(5, 0, 0), include.mean = TRUE, method = "ML")
# BIC values
bic_values <- c(
AR1 = Ar1$bic,
AR2 = Ar2$bic,
AR3 = Ar3$bic,
AR4 = Ar4$bic,
AR5 = Ar5$bic
)
# data frame
bic_results <- data.frame(
Model = names(bic_values),
BIC = bic_values
)
# Find the best model (lowest BIC)
best_model <- names(which.min(bic_values))
print("BIC Comparison Results:")
## [1] "BIC Comparison Results:"
print(bic_results)
## Model BIC
## AR1 AR1 3637.644
## AR2 AR2 3644.472
## AR3 AR3 3650.137
## AR4 AR4 3656.306
## AR5 AR5 3663.398
print(paste("The best model based on BIC is:", best_model))
## [1] "The best model based on BIC is: AR1"
# Extract AR(1) coefficients, standard errors, and t-statistics
ar1_coef <- coef(Ar1)
ar1_se <- sqrt(diag(Ar1$var.coef))
ar1_tstat <- ar1_coef / ar1_se
# AR(1) Forecast Evaluation
ar1f <- Arima(test_log_abs, model = Ar1)
fitted_ar1 <- fitted(ar1f)
# MSPE for log values
ar1_residuals_log <- test_log_abs - fitted_ar1
ar1_mspe_log <- sum(ar1_residuals_log^2) / length(ar1_residuals_log)
# MSPE for absolute values
ar1_residuals_exp <- exp(test_log_abs) - exp(fitted_ar1)
ar1_mspe_abs <- sum(ar1_residuals_exp^2) / length(ar1_residuals_exp)
# MEM(1,1) Model Estimation
g5 <- ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
distribution.model = "norm", fixed.pars = list(mu = 0)
)
mem11 <- ugarchfit(spec = g5, data = train_sqrt_abs)
# Extract MEM(1,1) coefficients, standard errors, and t-statistics
mem11_coef <- coef(mem11)
mem11_se <- mem11@fit$robust.se.coef
# MEM(1,1) Forecast Evaluation
g5f <- getspec(mem11)
setfixed(g5f) <- as.list(coef(mem11))
mem11_forecast <- ugarchforecast(g5f, n.ahead = 1, n.roll = 1040, data = test_sqrt_abs, out.sample = 1040)
mem_fitted <- mem11_forecast@forecast$sigmaFor[1, 1:1041]
# MSPE for absolute values
mem_residuals_abs <- test_sqrt_abs^2 - mem_fitted^2
mem_mspe_abs <- sum(mem_residuals_abs^2) / length(mem_residuals_abs)
# MSPE for log values
mem_residuals_log <- log(test_sqrt_abs^2) - log(mem_fitted^2)
mem_mspe_log <- sum(mem_residuals_log^2) / length(mem_residuals_log)
# Combine AR(1) Estimation Results
ar1_results <- data.frame(
Model = "AR(1)",
Parameter = names(ar1_coef),
Coefficient = ar1_coef,
SE = ar1_se,
T_Statistic = ar1_tstat
)
# Combine MEM(1,1) Estimation Results
mem11_results <- data.frame(
Model = rep("MEM(1,1)", length(mem11_coef)),
Parameter = names(mem11_coef),
Coefficient = mem11_coef,
SE = c(NA, mem11_se),
T_Statistic = mem11_coef / mem11_se
)
# Remove rows with NA
mem11_results <- mem11_results[complete.cases(mem11_results), ]
# Combine Estimation Results
estimation_results <- rbind(ar1_results, mem11_results)
# Combine Forecast Evaluation Results
forecast_evaluation <- data.frame(
Model = c("AR(1)", "MEM(1,1)"),
MSPE_abs = c(ar1_mspe_abs, mem_mspe_abs),
MSPE_log = c(ar1_mspe_log, mem_mspe_log)
)
print(estimation_results)
## Model Parameter Coefficient SE T_Statistic
## ar1 AR(1) ar1 -0.004643136 0.025834595 -1.797255e-01
## intercept AR(1) intercept -4.993732975 0.020755556 -2.405974e+02
## omega MEM(1,1) omega 0.000271959 0.002350849 1.626549e-03
## alpha1 MEM(1,1) alpha1 0.038632063 0.167200047 9.048045e-02
## beta1 MEM(1,1) beta1 0.930933696 0.426965855 3.959989e+02
print(forecast_evaluation)
## Model MSPE_abs MSPE_log
## 1 AR(1) 1.049185e-04 0.7940008
## 2 MEM(1,1) 7.454451e-05 0.7434122
Problem 5
# Step 1: Define binary indicator (ind) based on returns (ret)
ind <- ret
ind[ind > 0] <- 1
ind[!ind > 0] <- 0
# Step 2: Add lagged variables to data
dat_1$ind <- c(NaN, ind)
dat_1$ret <- c(NaN, ret)
dat_1$ret_l1 <- c(NaN, dat_1$ret[1:2541])
dat_1$ind_l1 <- c(NaN, dat_1$ind[1:2541])
# Step 3: Create train and test datasets
logit_data <- data.frame(
ind = dat_1$ind[3:nrow(dat_1)],
ret_l1 = dat_1$ret_l1[3:nrow(dat_1)],
ind_l1 = dat_1$ind_l1[3:nrow(dat_1)]
)
logit_train <- logit_data[1:1500, ]
logit_test <- logit_data[1501:nrow(logit_data), ]
# Step 4: Logistic regression models
logit_ret <- glm(ind ~ 1 + ret_l1, data = logit_train, family = "binomial") # Logit based on lagged returns
logit_ind <- glm(ind ~ 1 + ind_l1, data = logit_train, family = "binomial") # Logit based on lagged indicator
# Step 5: Create table for logit estimation results
logit_res <- data.frame(
Model = c("Lagged Returns (Logit)", "Lagged Indicator (Logit)"),
Intercept = c(coef(logit_ret)[1], coef(logit_ind)[1]),
RetL1 = c(coef(logit_ret)[2], NA),
IndL1 = c(NA, coef(logit_ind)[2]),
AIC = c(AIC(logit_ret), AIC(logit_ind))
)
print(logit_res)
## Model Intercept RetL1 IndL1 AIC
## ret_l1 Lagged Returns (Logit) 0.01565250 1.069839 NA 2083.298
## Lagged Indicator (Logit) -0.01075279 NA 0.05308716 2083.081
# Step 6: Make predictions and calculate accuracy
logit_test$ind_ret_pred <- predict(logit_ret, newdata = logit_test, type = "response")
logit_test$ind_ret_pred[logit_test$ind_ret_pred >= 0.5] <- 1
logit_test$ind_ret_pred[!logit_test$ind_ret_pred >= 0.5] <- 0
logit_test$ind_ind_pred <- predict(logit_ind, newdata = logit_test, type = "response")
logit_test$ind_ind_pred[logit_test$ind_ind_pred >= 0.5] <- 1
logit_test$ind_ind_pred[!logit_test$ind_ind_pred >= 0.5] <- 0
# Baseline models (always up and always down)
logit_test$up_pred <- rep(1, nrow(logit_test))
logit_test$down_pred <- rep(0, nrow(logit_test))
# Step 7: Calculate accuracy
logit_ret_acc <- sum(logit_test$ind_ret_pred == logit_test$ind) / nrow(logit_test)
logit_ind_acc <- sum(logit_test$ind_ind_pred == logit_test$ind) / nrow(logit_test)
up_acc <- sum(logit_test$up_pred == logit_test$ind) / nrow(logit_test)
down_acc <- sum(logit_test$down_pred == logit_test$ind) / nrow(logit_test)
# Step 8: Create forecast evaluation table
forecast_eval <- data.frame(
Model = c("Lagged Returns (Logit)", "Lagged Indicator (Logit)", "Always Up", "Always Down"),
Sign_criterion = c(logit_ret_acc, logit_ind_acc, up_acc, down_acc)
)
print(forecast_eval)
## Model Sign_criterion
## 1 Lagged Returns (Logit) 0.5500000
## 2 Lagged Indicator (Logit) 0.5365385
## 3 Always Up 0.5490385
## 4 Always Down 0.4509615