Factors that influence the volatility of stock prices

cvfitABSalpha0.2<-cv.biglasso(finalcut, yABS, alpha = 0.2)

elasticnetABSalpha0.2 <- biglasso(finalcut, yABS, alpha = 0.2 , lambda = cvfitABSalpha0.2$lambda.min)

cvfitABSalpha0.3<-cv.biglasso(finalcut, yABS, alpha = 0.3)

elasticnetABSalpha0.3 <- biglasso(finalcut, yABS, alpha = 0.3 , lambda = cvfitABSalpha0.3$lambda.min)

cvfitABSalpha0.4<-cv.biglasso(finalcut, yABS, alpha = 0.4)

elasticnetABSalpha0.4 <- biglasso(finalcut, yABS, alpha = 0.4 , lambda = cvfitABSalpha0.4$lambda.min)

cvfitABSalpha0<-cv.biglasso(finalcut, yABS, alpha = 0)

elasticnetABSalpha0 <- biglasso(finalcut, yABS, alpha = 0 , lambda = cvfitABSalpha0$lambda.min)

cvfitABSalpha0.6<-cv.biglasso(finalcut, yABS, alpha = 0.6)

elasticnetABSalpha0.6 <- biglasso(finalcut, yABS, alpha = 0.6 , lambda = cvfitABSalpha0.6$lambda.min)

cvfitABSalpha0.7<-cv.biglasso(finalcut, yABS, alpha = 0.7)

elasticnetABSalpha0.7 <- biglasso(finalcut, yABS, alpha = 0.7 , lambda = cvfitABSalpha0.7$lambda.min)

cvfitABSalpha0.8<-cv.biglasso(finalcut, yABS, alpha = 0.8)

elasticnetABSalpha0.8 <- biglasso(finalcut, yABS, alpha = 0.8 , lambda = cvfitABSalpha0.8$lambda.min)

cvfitABSalpha0.9<-cv.biglasso(finalcut, yABS, alpha = 0.9)

elasticnetABSalpha0.9 <- biglasso(finalcut, yABS, alpha = 0.9 , lambda = cvfitABSalpha0.9$lambda.min)

cvfitABSalpha1<-cv.biglasso(finalcut, yABS, alpha = 1)

elasticnetABSalpha1 <- biglasso(finalcut, yABS, alpha = 1 , lambda = cvfitABSalpha1$lambda.min)

ELASTICNETbetaABS=matrix(nrow=357, ncol=11)

ELASTICNETbetaABS[,1] <- elasticnetABSalpha0$beta[,1]

ELASTICNETbetaABS[,2] <- elasticnetABSalpha0.1$beta[,1]

ELASTICNETbetaABS[,3]<- elasticnetABSalpha0.2$beta[,1]

ELASTICNETbetaABS[,4] <- elasticnetABSalpha0.3$beta[,1]

ELASTICNETbetaABS[,5] <- elasticnetABSalpha0.4$beta[,1]

ELASTICNETbetaABS[,6] <- elasticnetABS$beta[,1]

ELASTICNETbetaABS[,7]<- elasticnetABSalpha0.6$beta[,1]

ELASTICNETbetaABS[,8] <- elasticnetABSalpha0.7$beta[,1]

ELASTICNETbetaABS[,9] <- elasticnetABSalpha0.8$beta[,1]

ELASTICNETbetaABS[,10] <- elasticnetABSalpha0.9$beta[,1]

ELASTICNETbetaABS[,11] <- elasticnetABSalpha1$beta[,1]

cvfitSQ <- cv.biglasso(finalcut,ySQ, eval.metric = c("default", "MAPE"), alpha = 0.5)

elasticnetSQ <- biglasso(finalcut, ySQ, penalty = c("enet"), alpha = 0.5, lambda = cvfitSQ$lambda.min)

cvfitSQalpha0.1<- cv.biglasso(finalcut, ySQ, eval.metric = c("default"), alpha = 0.1)

elasticnetSQalpha0.1 <- biglasso(finalcut, ySQ, alpha = 0.1 , lambda = cvfitSQalpha0.1$lambda.min)

cvfitSQalpha0.2<-cv.biglasso(finalcut, ySQ, alpha = 0.2)

elasticnetSQalpha0.2 <- biglasso(finalcut, ySQ, alpha = 0.2 , lambda = cvfitSQalpha0.2$lambda.min)

cvfitSQalpha0.3<-cv.biglasso(finalcut, ySQ, alpha = 0.3)

elasticnetSQalpha0.3 <- biglasso(finalcut, ySQ, alpha = 0.3 , lambda = cvfitSQalpha0.3$lambda.min)

cvfitSQalpha0.4<-cv.biglasso(finalcut, ySQ, alpha = 0.4)

elasticnetSQalpha0.4 <- biglasso(finalcut, ySQ, alpha = 0.4 , lambda = cvfitSQalpha0.4$lambda.min)

cvfitSQalpha0<-cv.biglasso(finalcut, ySQ, alpha = 0)

elasticnetSQalpha0 <- biglasso(finalcut, ySQ, alpha = 0 , lambda = cvfitSQalpha0$lambda.min)

cvfitSQalpha0.6<-cv.biglasso(finalcut, ySQ, alpha = 0.6)

elasticnetSQalpha0.6 <- biglasso(finalcut, ySQ, alpha = 0.6 , lambda = cvfitSQalpha0.6$lambda.min)

cvfitSQalpha0.7<-cv.biglasso(finalcut, ySQ, alpha = 0.7)

elasticnetSQalpha0.7 <- biglasso(finalcut, ySQ, alpha = 0.7 , lambda = cvfitSQalpha0.7$lambda.min)

cvfitSQalpha0.8<-cv.biglasso(finalcut, ySQ, alpha = 0.8)

elasticnetSQalpha0.8 <- biglasso(finalcut, ySQ, alpha = 0.8 , lambda = cvfitSQalpha0.8$lambda.min)

cvfitSQalpha0.9<-cv.biglasso(finalcut, ySQ, alpha = 0.9)

elasticnetSQalpha0.9 <- biglasso(finalcut, ySQ, alpha = 0.9 , lambda = cvfitSQalpha0.9$lambda.min)

cvfitSQalpha1<-cv.biglasso(finalcut, ySQ, alpha = 1)

elasticnetSQalpha1 <- biglasso(finalcut, ySQ, alpha = 1 , lambda = cvfitSQalpha1$lambda.min)

ELASTICNETbetaSQ=matrix(nrow=357, ncol=11)

ELASTICNETbetaSQ[,1] <- elasticnetSQalpha0$beta[,1]

ELASTICNETbetaSQ[,2] <- elasticnetSQalpha0.1$beta[,1]

ELASTICNETbetaSQ[,3]<- elasticnetSQalpha0.2$beta[,1]

ELASTICNETbetaSQ[,4] <- elasticnetSQalpha0.3$beta[,1]

ELASTICNETbetaSQ[,5] <- elasticnetSQalpha0.4$beta[,1]

ELASTICNETbetaSQ[,6] <- elasticnetSQ$beta[,1]

ELASTICNETbetaSQ[,7]<- elasticnetSQalpha0.6$beta[,1]

ELASTICNETbetaSQ[,8] <- elasticnetSQalpha0.7$beta[,1]

ELASTICNETbetaSQ[,9] <- elasticnetSQalpha0.8$beta[,1]

ELASTICNETbetaSQ[,10] <- elasticnetSQalpha0.9$beta[,1]

ELASTICNETbetaSQ[,11] <- elasticnetSQalpha1$beta[,1]

cvfitRB <- cv.biglasso(finalcut,yRB, eval.metric = c("default", "MAPE"), alpha = 0.5)

elasticnetRB <- biglasso(finalcut, yRB, penalty = c("enet"), alpha = 0.5, lambda = cvfitRB$lambda.min)

cvfitRBalpha0.1<-cv.biglasso(finalcut, yRB, alpha = 0.1)

elasticnetRBalpha0.1 <- biglasso(finalcut, yRB, alpha = 0.1 , lambda = cvfitRBalpha0.1$lambda.min)

cvfitABSalpha0.2<-cv.biglasso(finalcut, yABS, alpha = 0.2)

elasticnetRBalpha0.2 <- biglasso(finalcut, yRB, alpha = 0.2 , lambda = cvfitRBalpha0.2$lambda.min)

cvfitRBalpha0.3<-cv.biglasso(finalcut, yRB, alpha = 0.3)

elasticnetRBalpha0.3 <- biglasso(finalcut, yRB, alpha = 0.3 , lambda = cvfitRBalpha0.3$lambda.min)

cvfitRBalpha0.4<-cv.biglasso(finalcut, yRB, alpha = 0.4)

elasticnetRBalpha0.4 <- biglasso(finalcut, yRB, alpha = 0.4 , lambda = cvfitRBalpha0.4$lambda.min)

cvfitRBalpha0<-cv.biglasso(finalcut, yRB, alpha = 0)

elasticnetRBalpha0 <- biglasso(finalcut, yRB, alpha = 0 , lambda = cvfitRBalpha0$lambda.min)

cvfitRBalpha0.6<-cv.biglasso(finalcut, yRB, alpha = 0.6)

elasticnetRBalpha0.6 <- biglasso(finalcut, yRB, alpha = 0.6 , lambda = cvfitRBalpha0.6$lambda.min)

cvfitRBalpha0.7<-cv.biglasso(finalcut, yRB, alpha = 0.7)

elasticnetRBalpha0.7 <- biglasso(finalcut, yRB, alpha = 0.7 , lambda = cvfitRBalpha0.7$lambda.min)

cvfitRBalpha0.8<-cv.biglasso(finalcut, yRB, alpha = 0.8)

elasticnetRBalpha0.8 <- biglasso(finalcut, yRB, alpha = 0.8 , lambda = cvfitRBalpha0.8$lambda.min)

cvfitRBalpha0.9<-cv.biglasso(finalcut, yRB, alpha = 0.9)

elasticnetRBalpha0.9 <- biglasso(finalcut, yRB, alpha = 0.9 , lambda = cvfitRBalpha0.9$lambda.min)

cvfitRBalpha1<-cv.biglasso(finalcut, yRB, alpha = 1)

elasticnetRBalpha1 <- biglasso(finalcut, yRB, alpha = 1 , lambda = cvfitRBalpha1$lambda.min)

ELASTICNETbetaRB = matrix(nrow=357, ncol=11)

ELASTICNETbetaRB[,1] <- elasticnetRBalpha0$beta[,1]

ELASTICNETbetaRB[,2] <- elasticnetRBalpha0.1$beta[,1]

ELASTICNETbetaRB[,3]<- elasticnetRBalpha0.2$beta[,1]

ELASTICNETbetaRB[,4] <- elasticnetRBalpha0.3$beta[,1]

ELASTICNETbetaRB[,5] <- elasticnetRBalpha0.4$beta[,1]

ELASTICNETbetaRB[,6] <- elasticnetRB$beta[,1]

ELASTICNETbetaRB[,7]<- elasticnetRBalpha0.6$beta[,1]

ELASTICNETbetaRB[,8] <- elasticnetRBalpha0.7$beta[,1]

ELASTICNETbetaRB[,9] <- elasticnetRBalpha0.8$beta[,1]

ELASTICNETbetaRB[,10] <- elasticnetRBalpha0.9$beta[,1]

ELASTICNETbetaRB[,11] <- elasticnetRBalpha1$beta[,1]

grid <- seq( from = 0, to = 1, by = 0.1)

colnames(ELASTICNETbetaABS) <- grid

Intercept <- data.frame(filtered_names="INTERCEPT",filtered_class= 0)

names <- rbind(Intercept,names_f)

rownames(ELASTICNETbetaABS) <- names$filtered_names

colnames(ELASTICNETbetaSQ) <- grid

rownames(ELASTICNETbetaSQ) <- names$filtered_names

colnames(ELASTICNETbetaRB) <- grid

rownames(ELASTICNETbetaRB) <- names$filtered_names

colnames(ELASTICNETbetaVAR) <- grid

rownames(ELASTICNETbetaVAR) <- names$filtered_names

library(ddalpha)

library(caret)

library(lars)

library(tidyverse)

#larscvABS <- trainControl(finalcut, yABS , method = "lars")

print(lars)

#Linear Regression with Stepwise Selection (method = 'leapSeq')

#Linear Regression with Forward Selection (method = 'leapForward')

#Least Angle Regression (method = 'lars')For regression using package lars with tuning parameters:Fraction (fraction)

#Least Angle Regression (method = 'lars2') For regression using package lars with tuning parameters:Number of Steps (step)

ctrl <- trainControl(method = "cv", savePred=T, number=T)

#fold.ids <- createMultiFolds(yABS,k=10,times=3)

#fitControl <- trainControl(method = "repeatedcv", number = 10,repeats = 3,returnResamp = "final",index = fold.ids,summaryFunction = defaultSummary,selectionFunction = "oneSE")

yandxABC <- cbind(yABS, finalcutNOTBIG)

yandxABC <- as.data.frame(lapply(yandxABC, function(x) as.numeric(as.character(x))))

yandxSQ <- cbind(ySQ, finalcutNOTBIG)

yandxSQ <- as.data.frame(lapply(yandxSQ, function(x) as.numeric(as.character(x))))

yandxVAR <- cbind(yVARna, finalcutNOTBIGVAR)

yandxVAR <- as.data.frame(lapply(yandxVAR, function(x) as.numeric(as.character(x))))

yandxRB <- cbind(yRB, finalcutNOTBIG)

yandxRB <- as.data.frame(lapply(yandxRB, function(x) as.numeric(as.character(x))))

#ctrl <- trainControl(method = "cv", savePred=T, number=T)

library(xgboost)

set.seed(100) # For reproducibility

inTrainABS <- createDataPartition(y = yandxABC$weekly.returns, list = FALSE) #p = 0.8

trainingABS <- yandxABC[inTrainABS, ]

testingABS <- yandxABC[-inTrainABS,]

#for ABS

X_trainABS = xgb.DMatrix(as.matrix(trainingABS %>% select(-weekly.returns)))

y_trainABS = trainingABS$weekly.returns

X_testABS = xgb.DMatrix(as.matrix(testingABS %>% select(-weekly.returns)))

y_testABS = testingABS$weekly.returns

xgb_trcontrol = trainControl(method = "cv",number = 5,allowParallel = TRUE,verboseIter = FALSE,returnData = FALSE, initialWindow = 240,

horizon = 1,

fixedWindow = FALSE)

# by default parameters

xgbGrid <- expand.grid(nrounds = c(200), max_depth = c(15), colsample_bytree = (0.6),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelABS = train( X_trainABS, y_trainABS,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbTree")

#best valuesfor hyperparameters

xgb_modelABS$bestTune

xgbGridbest <- expand.grid(nrounds = 200, max_depth = 15, colsample_bytree = 0.6,eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelABSbest = train( X_trainABS, y_trainABS,trControl = xgb_trcontrol,tuneGrid = xgbGridbest, method = "xgbTree")

predictedABS = predict(xgb_modelABS, X_testABS) #or

#predictedABS = predict(xgb_modelbestABS, X_testABS)

residualsABS = y_testABS - predictedABS

RMSEabs = sqrt(mean(residualsABS^2))

cat('The root mean square error of the test data is ', round(RMSEabs,3),'\n')

#the root mean square error of the test data is 353.346

y_test_meanABS = mean(y_testABS)

# Calculate total sum of squares

tssABS = sum((y_testABS - y_test_meanABS)^2 )

# Calculate residual sum of squares

rssABS = sum(residualsABS^2)

# Calculate R-squared

rsqABS = 1 - (rssABS/tssABS)

cat('The R-square of the test data is ', round(rsqABS,3), '\n')

#The R-square of the test data is 0.129

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataABS = as.data.frame(cbind(predicted = predictedABS,

observed = y_testABS))

# Plot predictions vs test data

ggplot(my_dataABS,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#OR IF XGBLINEAR

xgbGrid <- expand.grid(nrounds = c(100,200), lambda=0,eta = 0.1, alpha=1)

xgbABSlin_model = train( X_trainABS, y_trainABS,trControl = xgb_trcontrol,tuneGrid = xgbGrid, label = TRUE,method = "xgbLinear")

xgbABSlin_model$bestTune

#xgbABSlinGridbest <- expand.grid(nrounds = 200, lambda = 0 , eta = 0.1, alpha =1 )

#xgb_modelABSlinbest = train( X_trainABS, y_trainABS,trControl = xgb_trcontrol,tuneGrid = xgbABSlinGridbest, method = "xgbLinear")

predictedABSlin = predict(xgbABSlin_model, X_testABS) #or

#predictedABSlin = predict(xgb_modelbestABSlinbest, X_testABS)

residualsABSlin = y_testABS - predictedABSlin

RMSEabslin = sqrt(mean(residualsABSlin^2))

cat('The root mean square error of the test data is ', round(RMSEabslin,3),'\n')

#the root mean square error of the test data is 148.203

y_test_meanABS = mean(y_testABS)

# Calculate total sum of squares

tssABS = sum((y_testABS - y_test_meanABS)^2 )

# Calculate residual sum of squares

rssABSlin = sum(residualsABSlin^2)

# Calculate R-squared

rsqABSlin = 1 - (rssABS/tssABS)

cat('The R-square of the test data is ', round(rsqABSlin,3),'\n')

#The R-square of the test data is 0.847

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataABSlin = as.data.frame(cbind(predicted = predictedABSlin,

observed = y_testABS))

# Plot predictions vs test data

ggplot(my_dataABSlin,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#for SQ

inTrainSQ <- createDataPartition(y = yandxSQ$weekly.returns, list = FALSE) #p = 0.8

trainingSQ <- yandxSQ[inTrainSQ, ]

testingSQ <- yandxSQ[-inTrainSQ,]

X_trainSQ = xgb.DMatrix(as.matrix(trainingSQ %>% select(-weekly.returns)))

y_trainSQ = trainingSQ$weekly.returns

X_testSQ = xgb.DMatrix(as.matrix(testingSQ %>% select(-weekly.returns)))

y_testSQ = testingSQ$weekly.returns

xgb_trcontrol = trainControl(method = "cv",number = 5,allowParallel = TRUE,verboseIter = FALSE,returnData = FALSE)

# by default parameters

xgbGrid <- expand.grid(nrounds = c(100,200), max_depth = c(10, 15, 20, 25), colsample_bytree = seq(0.5, 0.9, length.out = 5),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelSQ = train( X_trainSQ, y_trainSQ,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbTree", label = yandxABC$s)

#best valuesfor hyperparameters

xgb_modelSQ$bestTune #

xgbGridbest <- expand.grid(nrounds = c(200), max_depth = c(10), colsample_bytree = seq(0.8),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelSQbest = train( X_trainSQ, y_trainSQ,trControl = xgb_trcontrol,tuneGrid = xgbGridbest, method = "xgbTree")

predictedSQ = predict(xgb_modelSQ, X_testSQ) #or

#predictedSQ = predict(xgb_modelSQbest, X_testSQ)

residualsSQ = y_testSQ - predictedSQ

RMSEsq = sqrt(mean(residualsSQ^2))

cat('The root mean square error of the test data is ', round(RMSEsq,3),'\n')

#the root mean square error of the test data is 0.002

y_test_meanSQ = mean(y_testSQ)

# Calculate total sum of squares

tssSQ = sum((y_testSQ - y_test_meanSQ)^2 )

# Calculate residual sum of squares

rssSQ = sum(residualsSQ^2)

# Calculate R-squared

rsqSQ = 1 - (rssSQ/tssSQ)

cat('The R-square of the test data is ', round(rsqSQ,3), '\n')

#The R-square of the test data is 0.151

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataSQ = as.data.frame(cbind(predicted = predictedSQ,

observed = y_testSQ))

# Plot predictions vs test data

ggplot(my_dataSQ,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#OR IF XGBLINEAR

xgbGrid <- expand.grid(nrounds = c(100,200), lambda=0,eta = 0.1, alpha=1)

xgbSQlin_model = train( X_trainSQ, y_trainSQ,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbLinear")

xgbSQlin_model$bestTune

#xgbSQlinGridbest <- expand.grid(nrounds = 100, lambda =0 , eta = 0.1, alpha =1 )

#xgb_modelSQlinbest = train( X_trainABS, y_trainABS,trControl = xgb_trcontrol,tuneGrid = xgbABSlinGridbest, method = "xgbLinear")

predictedSQlin = predict(xgbSQlin_model, X_testSQ) #or

#predictedSQlin = predict(xgb_modelSQlinbest, X_testSQ)

residualsSQlin = y_testSQ - predictedSQlin

RMSEsqlin = sqrt(mean(residualsSQlin^2))

cat('The root mean square error of the test data is ', round(RMSEsqlin,3),'\n')

#the root mean square error of the test data is 0.003

y_test_meanSQ = mean(y_testSQ)

# Calculate total sum of squares

tssSQ = sum((y_testSQ - y_test_meanSQ)^2 )

# Calculate residual sum of squares

rssSQlin = sum(residualsSQlin^2)

# Calculate R-squared

rsqSQlin = 1 - (rssSQlin/tssSQ)

cat('The R-square of the test data is ', round(rsqSQlin,3), '\n')

#The R-square of the test data is

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataSQlin = as.data.frame(cbind(predicted = predictedSQlin,

observed = y_testSQ))

# Plot predictions vs test data

ggplot(my_dataSQlin,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#for RB

inTrainRB <- createDataPartition(y = yandxRB$SPY.High, list = FALSE) #p = 0.8

trainingRB <- yandxRB[inTrainRB, ]

testingRB <- yandxRB[-inTrainRB,]

X_trainRB = xgb.DMatrix(as.matrix(trainingRB %>% select(-SPY.High)))

y_trainRB = trainingRB$SPY.High

X_testRB = xgb.DMatrix(as.matrix(testingRB %>% select(-SPY.High)))

y_testRB = testingRB$SPY.High

xgb_trcontrol = trainControl(method = "cv",number = 5,allowParallel = TRUE,verboseIter = FALSE,returnData = FALSE)

# by default parameters

xgbGrid <- expand.grid(nrounds = c(100,200), max_depth = c(10, 15, 20, 25), colsample_bytree = seq(0.5, 0.9, length.out = 5),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelRB = train( X_trainRB, y_trainRB,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbTree")

#best valuesfor hyperparameters

xgb_modelRB$bestTune

xgbGridbest <- expand.grid(nrounds = c(100), max_depth = c(10), colsample_bytree = seq(0.7),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelRBbest = train( X_trainRB, y_trainRB,trControl = xgb_trcontrol,tuneGrid = xgbGridbest, method = "xgbTree")

predictedRB = predict(xgb_modelRB, X_testRB) #or

predictedRB = predict(xgb_modelRBbest, X_testRB)

residualsRB = y_testRB - predictedRB

RMSErb = sqrt(mean(residualsRB^2))

cat('The root mean square error of the test data is ', round(RMSErb,3),'\n')

#the root mean square error of the test data is 1.765

varImp(xgb_modelRB)

y_test_meanRB = mean(y_testRB)

# Calculate total sum of squares

tssRB = sum((y_testRB - y_test_meanRB)^2 )

# Calculate residual sum of squares

rssRB = sum(residualsRB^2)

# Calculate R-squared

rsqRB = 1 - (rssRB/tssRB)

cat('The R-square of the test data is ', round(rsqRB,3), '\n')

#The R-square of the test data is 0.476

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataRB = as.data.frame(cbind(predicted = predictedRB,

observed = y_testRB))

# Plot predictions vs test data

ggplot(my_dataRB,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#OR IF XGBLINEAR

xgbGrid <- expand.grid(nrounds = c(100,200), lambda=0,eta = 0.1, alpha=1)

xgbRBlin_model = train( X_trainRB, y_trainRB,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbLinear")

xgbRBlin_model$bestTune

xgbRBlinGridbest <- expand.grid(nrounds = 100 ,lambda = 0 , eta =0.1 , alpha = 1)

xgb_modelRBlinbest = train( X_trainRB, y_trainRB,trControl = xgb_trcontrol,tuneGrid = xgbRBlinGridbest, method = "xgbLinear")

predictedRBlin = predict(xgbRBlin_model, X_testRB) #or

predictedRBlin = predict(xgb_modelRBlinbest, X_testRB)

residualsRBlin = y_testRB - predictedRBlin

RMSErblin = sqrt(mean(residualsRBlin^2))

cat('The root mean square error of the test data is ', round(RMSErblin,3),'\n')

#the root mean square error of the test data is 1.776

y_test_meanRB = mean(y_testRB)

# Calculate total sum of squares

tssRB = sum((y_testRB - y_test_meanRB)^2 )

# Calculate residual sum of squares

rssRBlin = sum(residualsRBlin^2)

# Calculate R-squared

rsqRBlin = 1 - (rssRBlin/tssRB)

cat('The R-square of the test data is ', round(rsqRBlin,3),'\n')

#The R-square of the test data is 0.47

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataRBlin = as.data.frame(cbind(predicted = predictedRBlin,

observed = y_testRB))

# Plot predictions vs test data

ggplot(my_dataRBlin,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#for VAR

inTrainVAR <- createDataPartition(y = yandxVAR$daily.returns, list = FALSE) #p = 0.8

trainingVAR <- yandxVAR[inTrainVAR, ]

testingVAR <- yandxVAR[-inTrainVAR,]

X_trainVAR = xgb.DMatrix(as.matrix(trainingVAR %>% select(-daily.returns)))

y_trainVAR = trainingVAR$daily.returns

X_testVAR = xgb.DMatrix(as.matrix(testingVAR %>% select(-daily.returns)))

y_testVAR = testingVAR$daily.returns

xgb_trcontrol = trainControl(method = "cv",number = 5,allowParallel = TRUE,verboseIter = FALSE,returnData = FALSE)

# by default parameters

xgbGrid <- expand.grid(nrounds = c(100,200), max_depth = c(10, 15, 20, 25), colsample_bytree = seq(0.5, 0.9, length.out = 5),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelVAR = train( X_trainVAR, y_trainVAR,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbTree")

#best valuesfor hyperparameters

xgb_modelVAR$bestTune

xgbGridbest <- expand.grid(nrounds = c(200), max_depth = c(10), colsample_bytree = seq(0.5),eta = 0.1, gamma=0, min_child_weight = 1,subsample = 1)

xgb_modelVARbest = train( X_trainVAR, y_trainVAR,trControl = xgb_trcontrol,tuneGrid = xgbGridbest, method = "xgbTree")

predictedVAR = predict(xgb_modelVAR, X_testVAR) #or

predictedVAR = predict(xgb_modelVARbest, X_testVAR)

residualsVAR = y_testVAR - predictedVAR

RMSEvar = sqrt(mean(residualsVAR^2))

cat('The root mean square error of the test data is ', round(RMSEvar,3),'\n')

#the root mean square error of the test data is 0

y_test_meanVAR = mean(y_testVAR)

# Calculate total sum of squares

tssVAR = sum((y_testVAR - y_test_meanVAR)^2 )

# Calculate residual sum of squares

rssVAR = sum(residualsVAR^2)

# Calculate R-squared

rsqVAR = 1 - (rssVAR/tssVAR)

cat('The R-square of the test data is ', round(rsqVAR,3), '\n')

#The R-square of the test data is 0.142

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataVAR = as.data.frame(cbind(predicted = predictedVAR,

observed = y_testVAR))

# Plot predictions vs test data

ggplot(my_dataVAR,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

#OR IF XGBLINEAR

xgbGrid <- expand.grid(nrounds = c(100,200), lambda=0,eta = 0.1, alpha=1)

xgbVARlin_model = train( X_trainVAR, y_trainVAR,trControl = xgb_trcontrol,tuneGrid = xgbGrid, method = "xgbLinear")

xgbVARlin_model$bestTune

#xgbVARlinGridbest <- expand.grid(nrounds = 100 , lambda = 0, eta =0.1 , alpha = 1)

#xgb_modelVARlinbest = train( X_trainVAR, y_trainVAR,trControl = xgb_trcontrol,tuneGrid = xgbVARlinGridbest, method = "xgbLinear")

predictedVARlin = predict(xgbVARlin_model, X_testVAR) #or

#predictedVARlin = predict(xgb_modelVARlinbest, X_testVAR)

residualsVARlin = y_testVAR - predictedVARlin

RMSEvarlin = sqrt(mean(residualsVARlin^2))

cat('The root mean square error of the test data is ', round(RMSEvarlin,3),'\n')

#the root mean square error of the test data is 0.002

y_test_meanVAR = mean(y_testVAR)

# Calculate total sum of squares

tssVAR = sum((y_testVAR - y_test_meanVAR)^2 )

# Calculate residual sum of squares

rssVARlin = sum(residualsVARlin^2)

# Calculate R-squared

rsqVARlin = 1 - (rssVARlin/tssVAR)

cat('The R-square of the test data is ', round(rsqVARlin,3),'\n')

#The R-square of the test data is -16.321

#PLOT ACTUAL WITH PREDICTED

options(repr.plot.width=8, repr.plot.height=4)

my_dataVARlin = as.data.frame(cbind(predicted = predictedVARlin,

observed = y_testVAR))

# Plot predictions vs test data

ggplot(my_dataVARlin,aes(predicted, observed)) + geom_point(color = "darkred", alpha = 0.5) +

geom_smooth(method=lm)+ ggtitle('Linear Regression ') + ggtitle("Extreme Gradient Boosting: Prediction vs Test Data") +

xlab("Predecited Volatility ") + ylab("Observed Volatility") +

theme(plot.title = element_text(color="darkgreen",size=16,hjust = 0.5),

axis.text.y = element_text(size=12), axis.text.x = element_text(size=12,hjust=.5),

axis.title.x = element_text(size=14), axis.title.y = element_text(size=14))

varimp <- varImp(xgb_modelSQ, scale = FALSE)

plot(varimp, top = 20)

varImp(xgb_modelVAR, scale = FALSE)

varImp(xgb_modelABS, scale = FALSE)

varImp(xgb_modelRB, scale = FALSE)

varImp(xgb_modelABS, scale = FALSE)

library(corrplot)

correlations = cor(yandxABC)

corrplot(correlations, method="color")

library(fGarch)

garch1 <- garchFit(~ garch(1,1), data = SPYWEEKRETURN)

coef(garch)

predict(garch)

install.packages("rugarch")

library(rugarch)

gspec.ru <- ugarchspec(mean.model=list(armaOrder=c(0,0)), distribution="std")

gfit.ru <- ugarchfit(gspec.ru, SPYWEEKRETURN, n.ahead=1, n.roll = 660)

garch <- ugarchforecast(gfit.ru)

results1 <- garch@forecast$seriesFor[1,] + garch@forecast$sigmaFor[1,]

results2 <- garch@forecast$seriesFor[1,] - garch@forecast$sigmaFor[1,]

#one step ahead forecast

garch@forecast

gfit.ru@fit

sigmagarch <- as.data.frame(gfit.ru@fit[["sigma"]])

sigmagarch <- as.data.frame(gfit.ru@fit[["sigma"]])[660:1315, ]

#with predicted rolling window from 660 obs as training set start

predictedABS <- predictedABS[660:1315, ]

gfit.ru <- ugarchfit(gspec.ru, SPYWEEKRETURN, initialWindow = 660,horizon = 1315, fixedWindow = 1)

#DIEBOLD MARIANO TEST

library(forecast)

PROXYMODUL1 <- PROXYMODUL[660:1315]

PROXYRB1 <- PROXYRB[660:1315]

PROXYSQUARED1 <- PROXYSQUARED[660:1315]

PROXYVARIANCE1 <- PROXYVARIANCEna[657:1312]

ForecastErrorsXgbABS <- (PROXYMODUL1 - predictedABSlin)

ForecastErrorsGARCH <- (PROXYMODUL1 - sigmagarch)#dm test in r require row forecast errors

#MSE = loss

dm.test(ForecastErrorsGARCH, ForecastErrorsXgbABS, alternative = c("greater"))

ForecastErrorsXgbRB <- (PROXYRB1 - predictedRB)

#Diebold-Mariano Test

ForecastErrorsXgbSQ <- (PROXYSQUARED1 - predictedSQ)

ForecastErrorsXgVAR <- (PROXYVARIANCE1 - predictedVAR)

Diebold-Mariano Test

data: ForecastErrorsGARCHForecastErrorsXgb

DM = -4.3037, Forecast horizon = 1, Loss function power = 2, p-value = 1

alternative hypothesis: greater # model 2 is more accurate than 1

Diebold-Mariano Test

data: ForecastErrorsGARCHForecastErrorsXgb

DM = -4.3037, Forecast horizon = 1, Loss function power = 2, p-value =

1.936e-05 = 0

alternative hypothesis: two.sided

dm.test(ForecastErrorsGARCH, ForecastErrorsXgbSQ, alternative = c("greater"))

Diebold-Mariano Test

data: ForecastErrorsGARCHForecastErrorsXgbRB

DM = -4.0937, Forecast horizon = 1, Loss function power = 2, p-value = 1

alternative hypothesis: greater

plot(PROXYMODUL)

plot(PROXYVARIANCE)

plot(PROXYSQUARED)

plot(PROXYRB)

predictedABS2 <- predictedABS^2

normspec <- ugarchspec(mean.model=list(armaOrder=c(1,0)), distribution="norm")

Nmodel = ugarchfit(normspec,SPYWEEKRETURN)

Nfit<- as.data.frame(Nmodel@fit[["sigma"]])[660:1315, ]

Nfit2 <- (Nfit)^2

PROXYMODUL2 <- PROXYMODUL1^2

dm.test(( PROXYMODUL2- Nfit2), (PROXYMODUL2 - predictedABS2 ), alternative = c("greater"), h = 1, power = 2)

#Diebold-Mariano Test

data: (PROXYMODUL2 - Nfit2)(PROXYMODUL2 - predictedABS2)

DM = -2.6094, Forecast horizon = 1, Loss function power = 2, p-value =

0.009278

alternative hypothesis: two.sided



Reference

1. Amendola, A., & Storti, G. (2008). A GMM procedure for combining volatility forecasts. Computational Statistics & Data Analysis, 52(6), 3047-3060.

2. Christiansen, c., Schmeling, m., & Schrimpf, a. (2012). A comprehensive look at financial volatility prediction by economic variables. Journal of Applied Econometrics, 27(6), 956-977.

3. Chinco, A. M., Clark-Joseph, A. D., & Ye, M. (2017). Sparse signals in the cross-section of returns (No. w23933). National Bureau of Economic Research.

4. Farmer, L., Schmidt, L., & Timmermann, A. (2018). Pockets of predictability.

5. Ghysels, E., Santa-Clara, P., & Valkanov, R. (2006). Predicting volatility: getting the most out of return data sampled at different frequencies. Journal of Econometrics, 131(1-2), 59-95.

6. Jon Kleinberg Himabindu Lakkaraju Jure Leskovec Jens Ludwig Sendhil Mullainathan (2017). HUMAN DECISIONS AND MACHINE PREDICTIONS

7. Malakhov D.I. (2017) Sign, volatility and returns: differences and commonalities in predictability and stability of factors. Manuscript.

8. Malkiel, B. G. (1995). Returns from investing in equity mutual funds 1971 to 1991. The Journal of finance, 50(2), 549-572.

9. Rapach, D., & Zhou, G. (2013). Forecasting stock returns. In Handbook of economic forecasting (Vol. 2, pp. 328-383).

10. Tianqi Chen & Carlos Guestrin. (2017). XGBoost: A Scalable Tree Boosting System

11. Zou, H. (2006). The adaptive lasso and its oracle properties. J. Amer. Statist. Assoc. 101 1418-1429.

12. Zou, H. and HASTIE, T. (2005). Regularization and variable selection via the elastic net. J. R. Stat. Soc. Ser. B Stat. Methodol. 67 301-320.

Размещено на Allbest.ru