.Rhistory

forecast_len<- whole_len- burning
egarch_vx_roll<- ugarchroll( spec = egarch_vx,
data= ret,
n.ahead = 1,
forecast.length = forecast_len,
refit.every = 5)
load('egarch_model')
load('egarch_model')
load('egarch_model')
egarch_model$spec
egarch_model$spec
egarch_model$plot
rm(list=ls())
library(lubridate)
setwd('/Users/Eric/Documents/Vol_prediction')
SPXdata<- read.csv('SPX_rvol.csv')
rownames(SPXdata)<- ymd( SPXdata$DATE)
SPXdata$SPX2.rvol<- sqrt(SPXdata$SPX2.rv)
## preliminary analysis for SPX2.rv
acf(SPXdata$SPX2.rvol)
acf(SPXdata$SPX2.r)
library(tseries)
adf.test(SPXdata$SPX2.rvol)
kpss.test(SPXdata$SPX2.rvol, null = 'Level')
# ADF and KPSS test get both rejected. It is a flag of LRD.
plot(SPXdata$DATE, SPXdata$SPX2.rvol, type = 'n')
lines( SPXdata$DATE, SPXdata$SPX2.rvol)
## Fact:
## 1. SPX2.r shows week auto corr and SPX2.rv shows strong LRD
## 2. SPX2.r accepts the hypo of stationary
Box.test(SPXdata$SPX2.r, lag =1 , type= 'Ljung-Box')
shapiro.test(SPXdata$SPX2.r)
plot(density(SPXdata$SPX2.r))
library(moments)
kurtosis(SPXdata$SPX2.r)
## model SPX2.r by t dist
library(MASS)
t.pars<-fitdistr(SPXdata$SPX2.r, densfun = 't', start= list(m=0,s= 0.01 ,df= 1))
plot(density(SPXdata$SPX2.r), xlim= c(-.1,.1), ylim=c(-1, 55) )
par(new=TRUE)
curve( dt( (x- t.pars$estimate[1])/t.pars$estimate[2],
df= t.pars$estimate[3])/ t.pars$estimate[2],
from= -.1,
to= .1, xlim= c(-.1,.1),
ylim=c (-1, 55),
col= 'green')
library(forecast)
auto.arima(SPXdata$SPX2.r)
## auto.arima indicates ARMA(2,0) model for SPX2.r
## So apply standard ARMA(2,0)-GARCH(1,1) model for SPX2.r
library(rugarch)
sgarch<- ugarchspec(mean.model = list( armaOrder= c(2,0),
include.mean= TRUE),
distribution.model = 'std')
sgarch_fitted<- ugarchfit(sgarch, data =  SPXdata$SPX2.r)
whole_len<- length( SPXdata$SPX2.r)
burning_len<- 1000
forecast_len<- whole_len- burning_len
ret<- data.frame(SPX2.r= SPXdata$SPX2.r)
rownames(ret)<- rownames(SPXdata)
# rolling estimation and forecasting
sgarch_roll<- ugarchroll(spec = sgarch,
data= ret,
n.ahead = 1,
forecast.length =forecast_len,
refit.every = 5)
rm(list=ls())
library(lubridate)
setwd('/Users/Eric/Documents/Vol_prediction')
SPXdata<- read.csv('SPX_rvol.csv')
rownames(SPXdata)<- ymd( SPXdata$DATE)
SPXdata$SPX2.rvol<- sqrt(SPXdata$SPX2.rv)
## preliminary analysis for SPX2.rv
acf(SPXdata$SPX2.rvol)
acf(SPXdata$SPX2.r)
library(tseries)
adf.test(SPXdata$SPX2.rvol)
kpss.test(SPXdata$SPX2.rvol, null = 'Level')
# ADF and KPSS test get both rejected. It is a flag of LRD.
plot(SPXdata$DATE, SPXdata$SPX2.rvol, type = 'n')
lines( SPXdata$DATE, SPXdata$SPX2.rvol)
## Fact:
## 1. SPX2.r shows week auto corr and SPX2.rv shows strong LRD
## 2. SPX2.r accepts the hypo of stationary
Box.test(SPXdata$SPX2.r, lag =1 , type= 'Ljung-Box')
shapiro.test(SPXdata$SPX2.r)
plot(density(SPXdata$SPX2.r))
library(moments)
kurtosis(SPXdata$SPX2.r)
## model SPX2.r by t dist
library(MASS)
t.pars<-fitdistr(SPXdata$SPX2.r, densfun = 't', start= list(m=0,s= 0.01 ,df= 1))
plot(density(SPXdata$SPX2.r), xlim= c(-.1,.1), ylim=c(-1, 55) )
par(new=TRUE)
curve( dt( (x- t.pars$estimate[1])/t.pars$estimate[2],
df= t.pars$estimate[3])/ t.pars$estimate[2],
from= -.1,
to= .1, xlim= c(-.1,.1),
ylim=c (-1, 55),
col= 'green')
library(forecast)
auto.arima(SPXdata$SPX2.r)
## auto.arima indicates ARMA(2,0) model for SPX2.r
## So apply standard ARMA(2,0)-GARCH(1,1) model for SPX2.r
library(rugarch)
sgarch<- ugarchspec(mean.model = list( armaOrder= c(2,0),
include.mean= TRUE),
distribution.model = 'std')
sgarch_fitted<- ugarchfit(sgarch, data =  SPXdata$SPX2.r)
whole_len<- length( SPXdata$SPX2.r)
burning_len<- 1000
forecast_len<- whole_len- burning_len
ret<- data.frame(SPX2.r= SPXdata$SPX2.r)
rownames(ret)<- rownames(SPXdata)
# rolling estimation and forecasting
sgarch_roll<- ugarchroll(spec = sgarch,
data= ret,
n.ahead = 1,
forecast.length =forecast_len,
refit.every = 5)
sgarch_roll
rm(list=ls())
library(lubridate)
setwd('/Users/Eric/Documents/Vol_prediction')
SPXdata<- read.csv('SPX_rvol.csv')
rownames(SPXdata)<- ymd( SPXdata$DATE)
SPXdata$SPX2.rvol<- sqrt(SPXdata$SPX2.rv)
## preliminary analysis for SPX2.rv
acf(SPXdata$SPX2.rvol)
acf(SPXdata$SPX2.r)
library(tseries)
adf.test(SPXdata$SPX2.rvol)
kpss.test(SPXdata$SPX2.rvol, null = 'Level')
# ADF and KPSS test get both rejected. It is a flag of LRD.
plot(SPXdata$DATE, SPXdata$SPX2.rvol, type = 'n')
lines( SPXdata$DATE, SPXdata$SPX2.rvol)
## Fact:
## 1. SPX2.r shows week auto corr and SPX2.rv shows strong LRD
## 2. SPX2.r accepts the hypo of stationary
Box.test(SPXdata$SPX2.r, lag =1 , type= 'Ljung-Box')
shapiro.test(SPXdata$SPX2.r)
plot(density(SPXdata$SPX2.r))
library(moments)
kurtosis(SPXdata$SPX2.r)
## model SPX2.r by t dist
library(MASS)
t.pars<-fitdistr(SPXdata$SPX2.r, densfun = 't', start= list(m=0,s= 0.01 ,df= 1))
plot(density(SPXdata$SPX2.r), xlim= c(-.1,.1), ylim=c(-1, 55) )
par(new=TRUE)
curve( dt( (x- t.pars$estimate[1])/t.pars$estimate[2],
df= t.pars$estimate[3])/ t.pars$estimate[2],
from= -.1,
to= .1, xlim= c(-.1,.1),
ylim=c (-1, 55),
col= 'green')
library(forecast)
auto.arima(SPXdata$SPX2.r)
## auto.arima indicates ARMA(2,0) model for SPX2.r
## So apply standard ARMA(2,0)-GARCH(1,1) model for SPX2.r
library(rugarch)
sgarch<- ugarchspec(mean.model = list( armaOrder= c(2,0),
include.mean= TRUE),
distribution.model = 'std')
sgarch_fitted<- ugarchfit(sgarch, data =  SPXdata$SPX2.r)
whole_len<- length( SPXdata$SPX2.r)
burning_len<- 1000
forecast_len<- whole_len- burning_len
ret<- data.frame(SPX2.r= SPXdata$SPX2.r)
rownames(ret)<- rownames(SPXdata)
# rolling estimation and forecasting
sgarch_roll<- ugarchroll(spec = sgarch,
data= ret,
n.ahead = 1,
forecast.length =forecast_len,
refit.every = 5)
## plot the predicted vol and realized vol(5-min estimator by OxManLab)
library(ggplot2)
library(reshape2)
x<- tail( ymd( SPXdata$DATE), forecast_len)
realized_vol<- sqrt(tail( SPXdata$SPX2.rv, forecast_len))
sgarch.predicted_vol<- sgarch_roll@forecast$density[,'Sigma']
tmp_df<- data.frame(x, realized_vol, sgarch.predicted_vol)
sgarch.g<- ggplot(melt(tmp_df, id.var= 'x'), aes(x=x, y= value))+
geom_line(aes(colour= variable, group= variable))+
scale_color_manual(values = c('grey', 'red'))+
ylab('daily volatility')+
xlab('date index')+
theme(legend.title= element_blank())+
ggtitle('ARMA(2,0)-GARCH(1,1) vol prediction')
jpeg('ARMAGARCH.jpeg')
sgarch.g
dev.off()
## compute the squared error
sgarch.MSE<- mean(( realized_vol- sgarch.predicted_vol)^2)
summary( (realized_vol-sgarch.predicted_vol)^2)
sgarch.MSE
cor( realized_vol, sgarch.predicted_vol)
summary( lm( realized_vol~ sgarch.predicted_vol))
## save sgarch model
sgarch_model<- list()
sgarch_model$spec<- sgarch
sgarch_model$roll<- sgarch_roll
sgarch_model$plot<- sgarch.g
sgarch_model$MSE<- sgarch.MSE
sgarch_model$roll.pred<- tmp_df
save(sgarch_model, file = 'sgarch_model')
## EGARCH estimation
egarch<- ugarchspec(variance.model = list( model='eGARCH'),
mean.model = list( armaOrder= c(2,0)),
distribution.model = 'std')
egarch_fitted<- ugarchfit(egarch, data= SPXdata$SPX2.r)
egarch_roll<- ugarchroll( spec = egarch,
data= ret,
n.ahead = 1,
forecast.length = forecast_len,
refit.every = 5)
x<- tail( ymd( SPXdata$DATE), forecast_len)
realized_vol<- sqrt(tail( SPXdata$SPX2.rv, forecast_len))
egarch.predicted_vol<- egarch_roll@forecast$density[,'Sigma']
tmp_df<- data.frame(x, realized_vol, egarch.predicted_vol)
egarch.g<- ggplot(melt(tmp_df, id.var= 'x'), aes(x=x, y= value))+
geom_line(aes(colour= variable, group= variable))+
scale_color_manual(values = c('grey', 'blue'))+
ylab('daily volatility')+
xlab('date index')+
theme(legend.title= element_blank())+
ggtitle('ARMA(2,0)-EGARCH(1,1) vol prediction')
jpeg('ARMAEGARCH.jpeg')
egarch.g
dev.off()
egarch.MSE<- mean(( realized_vol- egarch.predicted_vol)^2)
summary( (realized_vol-egarch.predicted_vol)^2)
egarch.MSE
cor( realized_vol, egarch.predicted_vol)
summary( lm( realized_vol~ egarch.predicted_vol))
egarch_model<- list()
egarch_model$spec<- egarch
egarch_model$roll<- egarch_roll
egarch_model$plot<- egarch.g
egarch_model$MSE<- egarch.MSE
egarch_model$roll.pred<- tmp_df
save(egarch_model, file= 'egarch_model')
load('egarch_model')
egarch_model$plot
cor( egarch_model$roll.pred$realized_vol, egarch_model$roll.pred$egarch.predicted_vol, method = 'spearman')
cor( egarch_model$roll.pred$realized_vol, egarch_model$roll.pred$egarch.predicted_vol, method = 'spearman')
cor( egarch_model$roll.pred$realized_vol, egarch_model$roll.pred$egarch.predicted_vol, method = 'spearman')
cor( egarch_model$roll.pred$realized_vol, egarch_model$roll.pred$egarch.predicted_vol, method = 'spearman')
summary(egarch_model$roll.pred$realized_vol-
egarch_model$roll.pred$egarch.predicted_vol)
rgarch_model$plot
rgarch_model$plot
load('rgarch_model')
rgarch_model$spec
rgarch_model$plot
rm(list=ls())
setwd('/Users/Eric/Desktop/Vol_prediction')
library(lubridate)
setwd('/Users/Eric/Desktop/Vol_prediction')
SPXdata<- read.csv('SPX_rvol.csv')
rownames(SPXdata)<- ymd( SPXdata$DATE)
SPXdata$SPX2.rvol<- sqrt(SPXdata$SPX2.rv)
library(rugarch)
library(xts)
rgarch.model<- ugarchspec(mean.model = list(armaOrder= c(2,0)),
variance.model = list(model= 'realGARCH',
garchOrder= c(2,1)))
setbounds(rgarch.model)<- list(alpha2=c (-1,1))
SPXdata.xts<- SPXdata
SPXdata.xts$DATE<- NULL
SPXdata.xts<- as.xts(SPXdata.xts)
rgarch.fit<- ugarchfit(spec = rgarch.model,
data= SPXdata.xts$SPX2.r,
solver= 'hybrid',
realizedVol= SPXdata.xts$SPX2.rvol)
whole_len<- dim(SPXdata.xts)[1]
burning<- 1000
forecast_len<- whole_len- burning
rgarch.roll<- ugarchroll(spec = rgarch.model,
data= SPXdata.xts$SPX2.r,
n.ahead = 1,
forecast.length = forecast_len,
refit.every = 5,
solver= 'hybrid',
realizedVol= SPXdata.xts$SPX2.rvol
)
library(ggplot2)
library(reshape2)
tmp_df<- data.frame(x= tail( ymd(index(SPXdata.xts)) ,forecast_len),
realized_vol= tail( SPXdata$SPX2.rvol, forecast_len),
rgarch.prediction_vol= rgarch.roll@forecast$density$RVolForecast)
rgarch.g<- ggplot( data= melt( tmp_df, id.var= 'x'), aes( x=x, y= value))+
geom_line(aes(colour= variable, group= variable))+
scale_colour_manual(values= c( 'grey', 'orange'))+
ylab('daily volatility')+
xlab( 'date')+
theme( legend.title= element_blank())+
ggtitle( 'realizedGARCH(ARMA(2,0)-rGARCH(2,1)) rolling prediction')
jpeg( 'rGARCH.jpeg')
rgarch.g
dev.off()
rgarch.MSE<- mean( (tmp_df$realized_vol- tmp_df$rgarch.prediction_vol)^2)
summary( (tmp_df$realized_vol- tmp_df$rgarch.prediction_vol)^2)
rgarch.MSE
cor( tmp_df$realized_vol, tmp_df$rgarch.prediction_vol)
summary( lm( tmp_df$realized_vol~ tmp_df$rgarch.prediction_vol))
rgarch_model<- list()
rgarch_model$spec<- rgarch.model
rgarch_model$roll<- rgarch.roll
rgarch_model$plot<- rgarch.g
rgarch_model$MSE<- rgarch.MSE
rgarch_model$roll.pred<- tmp_df
save(rgarch_model, file = 'rgarch_model')
rgarch_model$plot
rm(list=ls())
setwd('/Users/Eric/Desktop/Vol_prediction')
SPXdata<- read.csv('SPX_rvol.csv')
library(lubridate)
rownames(SPXdata)<- ymd(SPXdata$DATE)
SPXdata$SPX2.rvol<- sqrt(SPXdata$SPX2.rv)
acf( SPXdata$SPX2.rvol, 1000)
library(tseries)
library(pracma)
hurstexp(SPXdata$SPX2.rvol)
adf.test(SPXdata$SPX2.rvol)
kpss.test(SPXdata$SPX2.rvol)
## acf of SPX2.rv decays slowly. It indicates LRD.
## Model SPX2.rv by fractional ARIMA
library(rugarch)
## armaOrder should be (0,0)
# pq<- c()
# aic<- c()
# for (p in 0:1){
#   for (q in 0:1){
#     arfima.model<- ugarchspec(mean.model = list(armaOrder= c(p, q),
#                                                 arfima=TRUE))
#     arfima.fitted<- ugarchfit( spec= arfima.model,
#                                data= SPXdata$SPX2.rv)
#     pq<- c(pq, paste(c(p,q), collapse = ','))
#     aic<- c(aic, infocriteria(arfima.fitted)[1])
#   }
# }
# tmp.df<- data.frame(pq= pq, aic= aic)
# tmp.df
rVol<- data.frame(rvol= SPXdata$SPX2.rvol)
rownames(rVol)<- ymd(SPXdata$DATE)
arfima_egarch.model<- ugarchspec(mean.model = list(armaOrder= c(0,0),
arfima=TRUE),
variance.model = list(model='eGARCH'))
arfima_egarch.fitted<- ugarchfit(spec = arfima_egarch.model,
data= rVol)
arfima_egarch.fitted
library(ggplot2)
library(lubridate)
library(reshape2)
tmp_df<- data.frame(x= rownames(rVol),
realized= rVol$rvol,
arfima_egarch.fitted= arfima_egarch.fitted@fit$fitted.values )
# g<- ggplot(melt(tmp_df, id.vars = 'x'), aes(x=x, y= value))+
#   geom_line(aes(colour= variable, group= variable))+
#   scale_color_manual(values= c('grey', 'red'))+
#   ylab('daily volatility')+
#   xlab('date index')+
#   theme(legend.title= element_blank())
# g
cor( tmp_df$realized, tmp_df$arfima_egarch.fitted)
summary( lm( tmp_df$realized~ tmp_df$arfima_egarch.fitted))
##residual analysis
resi<- arfima_egarch.fitted@fit$residuals
plot(resi, type= 'l')
acf(resi, 30)
adf.test(resi)
kpss.test(resi)
library(lmtest)
## rolling prediction
whole_len<- dim(rVol)[1]
burning_len<- 1000
forecast_len<- whole_len- burning_len
arfima_egarch.roll<- ugarchroll(arfima_egarch.model,
data = rVol,
n.ahead = 1,
forecast.length = forecast_len,
refit.every = 5)
library(ggplot2)
library(reshape2)
tmp_df<- data.frame(x= tail( ymd( SPXdata$DATE), forecast_len),
realized_vol= tail(rVol$rvol, forecast_len),
arfima_egarch.predicted_vol= arfima_egarch.roll@forecast$density[,'Mu'])
arfima_egarch.g<- ggplot(melt(tmp_df, id.var= 'x'), aes(x=x, y= value))+
geom_line(aes(colour= variable, group= variable))+
scale_color_manual(values = c('grey', 'purple'))+
ylab('daily volatility')+
xlab('date index')+
theme(legend.title= element_blank())+
ggtitle('Realized Vol ARFIMA(0,d,0)_EGARCH(1,1) rolling prediction')
jpeg('rVolARFIMAEGARCH.jpeg')
arfima_egarch.g
dev.off()
## compute the squared error
arfima_egarch.MSE<- mean(( tmp_df$realized_vol- tmp_df$arfima_egarch.predicted_vol)^2)
summary( (tmp_df$realized_vol-tmp_df$arfima_egarch.predicted_vol)^2)
arfima_egarch.MSE
cor( tmp_df$realized_vol, tmp_df$arfima_egarch.predicted_vol)
summary( lm( tmp_df$realized_vol~ tmp_df$arfima_egarch.predicted_vol))
arfima_egarch_model<- list()
arfima_egarch_model$spec<- arfima_egarch.model
arfima_egarch_model$roll<- arfima_egarch.roll
arfima_egarch_model$plot<- arfima_egarch.g
arfima_egarch_model$MSE<- arfima_egarch.MSE
arfima_egarch_model$roll.pred<- tmp_df
save(arfima_egarch_model, file='arfima_egarch_model')
arfima_egarch_model$plot
arfima_egarch_model$plot
rm(list=ls())
setwd('/Users/Eric/Desktop/Vol_prediction')
load('egarch_model')
load('arfima_egarch_model')
load('rgarch_model')
library(lubridate)
library(xts)
library(zoo)
SPXdata<- read.csv('SPX_rvol.csv')
SPXdata$SPX2.rvol<- sqrt( SPXdata$SPX2.rv)
rownames(SPXdata)<- ymd( SPXdata$DATE)
## VX_lag1 contains the information of implied vol with 1 day lag
VX<- read.csv('VX.csv')
VX$DATE<- ymd(VX$DATE)
rownames(VX)<- ymd(VX$DATE)
VX$DATE<- NULL
VX<- as.xts(VX)
VX_lag1<- lag( VX, 1, na.pad = TRUE)
VX_lag1<- as.data.frame(VX_lag1)
VX_lag1$DATE<- ymd(rownames(VX_lag1))
VX_lag1<- VX_lag1[ complete.cases(VX_lag1), ]
## forecasting is the time series 1 day ahead forecasting
forecasting<- data.frame(DATE=ymd(egarch_model$roll.pred$x),
real= egarch_model$roll.pred$realized_vol,
egarchPred= egarch_model$roll.pred$egarch.predicted_vol,
rgarchPred= rgarch_model$roll.pred$rgarch.prediction_vol,
arfima_egarchPred= arfima_egarch_model$roll.pred$arfima_egarch.predicted_vol)
## combine the VX_lag1 and forecasting
fore_vx<- merge( VX_lag1, forecasting, by='DATE')
fore_vx<- fore_vx[complete.cases(fore_vx),]
rownames(fore_vx)<- ymd(fore_vx$DATE)
fore_vx$DATE<- NULL
fore_vx<- as.xts( fore_vx)
tmp<- list()
tmp$rvol<- fore_vx$real
fore_vx$real<- NULL
tmp$x<- fore_vx
fore_vx<- tmp
## rf ensemble to get an average forecasting ( not successful)
library(randomForest)
# # cv<- rfcv(trainx = fore_vx$x, trainy = fore_vx$rvol, cv.fold = 5)
#
# ensemble_rf<- randomForest(x = fore_vx$x ,
#                            y= fore_vx$rvol,
#                            ntree= 1000,
#                            importance= TRUE)
#
# plot(ensemble_rf)
# varImpPlot(ensemble_rf)
# var_selection<- c( 'rgarchPred',
#                    'arfima_egarchPred',
#                    'VIX.Close',
#                    'egarchPred',
#                    'VIX.Open',
#                    'VX.C5_2',
#                    'VX.C2_1',
#                    'VX.C1_0',
#                    'VX4.OpenInt_lag1',
#                    'VX1.Settle')
# tmp<- fore_vx$x[,  var_selection]
#
# ensemble_rf<- randomForest(x = tmp ,
#                            y= fore_vx$rvol,
#                            ntree= 1000,
#                            importance= TRUE)
#
## rf ensemble the time series forecasting to get an average
rownames(forecasting)<- forecasting$DATE
forecasting$DATE<- NULL
ensemble_rf<- randomForest(real~.,
data = forecasting,
ntree= 1000,
mtry= 2,
importance= TRUE)
ensemble_rf
plot(ensemble_rf)
ensemble_rf$mse
ensemble_rf$rsq
library(ggplot2)
library(reshape2)
tmp<- data.frame( x= ymd( egarch_model$roll.pred$x),
reallized_vol= forecasting$real,
rf.predicted_vol= ensemble_rf$predicted)
rf.g<- ggplot(melt(data =  tmp, id.var='x'), aes( x=x, y= value))+
geom_line(aes(colour= variable, group= variable))+
scale_color_manual(values = c('grey', 'dark green'))+
ylab('daily volatility')+
xlab('date index')+
theme(legend.title= element_blank())+
ggtitle('RandomForest ensemble prediction')
jpeg('RF_ENSEMBLE.jpeg')
rf.g
dev.off()
rf<- list()
rf$model<- ensemble_rf
rf$plot<- rf.g
rf$MSE<- mean( ensemble_rf$mse)
rf$roll.pred<- tmp
save(rf, file='rf')