#读入数据
library(tseries)
p = read.csv('/Users/lisiqi/Desktop/hs.csv')
p <- p$cls
y=diff(log(p))*100

#描述统计
summary(y)
sd(y)
library(moments)
skewness(y)
kurtosis(y)
jarque.bera.test(y)
adf.test(y)
t.test(y)

#时序图、ACF、PACF
plot(y,type="l",main="Time Series Plot of HS300.return",xlab="time",ylab="log return")
par(mfrow=c(2,1))
acf(y,main="ACF")
pacf(y,main="PACF")

#ARCH效应检验
mod_ar <- arima(y,order=c(3,0,0))
summary(mod_ar)
y <- mod_ar$residuals
Box.test(y^2,lag=18)

## We multiply returns by 100 and de-mean them
library(fGarch)
mod11 <- garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="norm",trace=F)
mod12 <- garchFit(~ garch(1,2), data = y,include.mean=FALSE,cond.dist="norm",trace=F)
mod13 <- garchFit(~ garch(2,1), data = y,include.mean=FALSE,cond.dist="norm",trace=F)
mod21 <- garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="std",trace=F)
mod22 <- garchFit(~ garch(1,2), data = y,include.mean=FALSE,cond.dist="std",trace=F)
mod23 <- garchFit(~ garch(2,1), data = y,include.mean=FALSE,cond.dist="std",trace=F)
mod31 <- garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="ged",trace=F)
mod32 <- garchFit(~ garch(1,2), data = y,include.mean=FALSE,cond.dist="ged",trace=F)
mod33 <- garchFit(~ garch(2,1), data = y,include.mean=FALSE,cond.dist="ged",trace=F)
summary(mod33)
mod1 <- garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="norm",trace=F)
mod2 <- garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="std",trace=F)
mod3 <- garchFit(~ garch(1,1), data = y,include.mean=FALSE,cond.dist="ged",trace=F)
summary(mod1)
summary(mod2)
summary(mod3)

#模型检验
r1 <- mod1@residuals
r2 <- mod2@residuals
r3 <- mod3@residuals
Box.test(r1,18)
Box.test(r2,18)
Box.test(r3,18)

#GARCH(1,1)模型预测VaR
library(tseries)
library(forecast)
library(fGarch)
library(fracdiff)                   
library(rugarch)
require(xts)
library(parallel)

p = read.csv('/Users/lisiqi/Desktop/hs.csv')
p <- p$cls
y=diff(log(p))*100
spec.garch1 = ugarchspec(mean.model=list(armaOrder=c(4,0),include.mean=TRUE),variance.model=list(model ='fGARCH',garchOrder=c(1,1),submodel='NAGARCH'),distribution.model="norm")
garchfit1 = ugarchfit(spec.garch1, y)
nagarch.norm1 = ugarchroll(spec.garch1,data=y,forecast.length=300,forerefit.window="recursive",window.size=1000,solver="hybrid",calculate.VaR= TRUE,VaR.alpha=c(0.01,0.05),keep.coef=TRUE)
VaR_g1 <- as.data.frame(nagarch.norm1,which="VaR")
var_95_1 <- VaR_g1$`alpha(5%)`
var_99_1 <- VaR_g1$`alpha(1%)`
mean(var_95_1)
mean(var_99_1)
sigma1 <- sigma(garchfit1)

spec.garch2 = ugarchspec(mean.model=list(armaOrder=c(4,0),include.mean=TRUE),variance.model=list(model ='fGARCH',garchOrder=c(1,1),submodel='NAGARCH'),distribution.model="std")
garchfit2 = ugarchfit(spec.garch2, y)
nagarch.norm2 = ugarchroll(spec.garch2,data=y,forecast.length=300,forerefit.window="recursive",window.size=1000,solver="hybrid",calculate.VaR= TRUE,VaR.alpha=c(0.01,0.05),keep.coef=TRUE)
VaR_g2 <- as.data.frame(nagarch.norm2,which="VaR")
var_95_2 <- VaR_g2$`alpha(5%)`
var_99_2 <- VaR_g2$`alpha(1%)`
mean(var_95_2)
mean(var_99_2)
sigma2 <- sigma(garchfit2)

spec.garch3 = ugarchspec(mean.model=list(armaOrder=c(4,0),include.mean=TRUE),variance.model=list(model ='fGARCH',garchOrder=c(1,1),submodel='NAGARCH'),distribution.model="ged")
garchfit3 = ugarchfit(spec.garch3, y)
nagarch.norm3 = ugarchroll(spec.garch3,data=y,forecast.length=300,forerefit.window="recursive",window.size=1000,solver="hybrid",calculate.VaR= TRUE,VaR.alpha=c(0.01,0.05),keep.coef=TRUE)
VaR_g3 <- as.data.frame(nagarch.norm3,which="VaR")
var_95_3 <- VaR_g3$`alpha(5%)`
var_99_3 <- VaR_g3$`alpha(1%)`
mean(var_95_3)
mean(var_99_3)
sigma3 <- sigma(garchfit3)

#分位数回归预测VaR
#分位数回归
n <- length(y)
data1 <- data.frame(y,sigma1,sigma1_2=sigma1^2)
data2 <- data.frame(y,sigma2,sigma2_2=sigma2^2)
data3 <- data.frame(y,sigma3,sigma3_2=sigma3^2)
#QR_GARCH(1,1)+norm
library(quantreg)
qr11 <- rq(y~sigma1+sigma1_2,data=data1,tau=0.05)
qr12 <- rq(y~sigma1+sigma1_2,data=data1,tau=0.01)
summary(qr11)
summary(qr12)

var_95_qr1 <- as.vector(predict(qr11,data1)[(n-299):n])
var_99_qr1 <- as.vector(predict(qr12,data1)[(n-299):n])
mean(var_95_qr1)
mean(var_99_qr1)

#QR_GARCH(1,1)+std
qr21 <- rq(y~sigma2+sigma2_2,data=data2,tau=0.05)
qr22 <- rq(y~sigma2+sigma2_2,data=data2,tau=0.01)
summary(qr21)
summary(qr22)

var_95_qr2 <- as.vector(predict(qr21,data2)[(n-299):n])
var_99_qr2 <- as.vector(predict(qr22,data2)[(n-299):n])
mean(var_95_qr2)
mean(var_99_qr2)

#QR_GARCH(1,1)+ged
qr31 <- rq(y~sigma3+sigma3_2,data=data3,tau=0.05)
qr32 <- rq(y~sigma3+sigma3_2,data=data3,tau=0.01)
summary(qr31)
summary(qr32)

var_95_qr3 <- as.vector(predict(qr31,data3)[(n-299):n])
var_99_qr3 <- as.vector(predict(qr32,data3)[(n-299):n])
mean(var_95_qr3)
mean(var_99_qr3)

#模型比较
#GARCH(1,1)+norm+95%VaR
par(mfrow=c(1,1))
plot(y[(n-299):n],xlab="time",ylab="log return",main="normal distribution 95% VaR")
lines(var_95_1,col="blue")
lines(var_95_qr1,col="red")
legend("topleft",c("GARCH","QR_GARCH"),col=c("blue","red"),lty=c(1,1))
sum(y[(n-299):n]<var_95_1)
sum(y[(n-299):n]<var_95_qr1)

#GARCH(1,1)+norm+99%VaR
plot(y[(n-299):n],xlab="time",ylab="log return",main="normal distribution 99% VaR")
lines(var_99_1,col="blue")
lines(var_99_qr1,col="red")
legend("topleft",c("GARCH","QR_GARCH"),col=c("blue","red"),lty=c(1,1))
sum(y[(n-299):n]<var_99_1)
sum(y[(n-299):n]<var_99_qr1)

#GARCH(1,1)+std+95%VaR
plot(y[(n-299):n],xlab="time",ylab="log return",main="student t distribution 95% VaR")
lines(var_95_2,col="blue")
lines(var_95_qr2,col="red")
legend("topleft",c("GARCH","QR_GARCH"),col=c("blue","red"),lty=c(1,1))
sum(y[(n-299):n]<var_95_2)
sum(y[(n-299):n]<var_95_qr2)

#GARCH(1,1)+std+99%VaR
plot(y[(n-299):n],xlab="time",ylab="log return",main="student t distribution 99% VaR")
lines(var_99_2,col="blue")
lines(var_99_qr2,col="red")
legend("topleft",c("GARCH","QR_GARCH"),col=c("blue","red"),lty=c(1,1))
sum(y[(n-299):n]<var_99_2)
sum(y[(n-299):n]<var_99_qr2)

#GARCH(1,1)+ged+95%VaR
plot(y[(n-299):n],xlab="time",ylab="log return",main="GED distribution 95% VaR")
lines(var_95_3,col="blue")
lines(var_95_qr3,col="red")
legend("topleft",c("GARCH","QR_GARCH"),col=c("blue","red"),lty=c(1,1))
sum(y[(n-299):n]<var_95_3)
sum(y[(n-299):n]<var_95_qr3)

#GARCH(1,1)+ged+99%VaR
plot(y[(n-299):n],xlab="time",ylab="log return",main="GED distribution 99% VaR")
lines(var_99_3,col="blue")
lines(var_99_qr3,col="red")
legend("topleft",c("GARCH","QR_GARCH"),col=c("blue","red"),lty=c(1,1))
sum(y[(n-299):n]<var_99_3)
sum(y[(n-299):n]<var_99_qr3)

#backtesting
m <- c(18,23,20,17,16,16)
n <- 300
LR1 <- -2*log(((1-0.05)^(n-m))*(0.05^m))+2*log(((1-m/n)^(n-m))*((m/n)^m))  #95%
LR1
1-pchisq(LR1,1)
m <- c(7,6,6,2,3,2)
LR2 <- -2*log(((1-0.01)^(n-m))*(0.01^m))+2*log(((1-m/n)^(n-m))*((m/n)^m))  #99%
LR2
1-pchisq(LR2,1)

#两市场比较
result_hs <- cbind(var_95_qr1,var_95_qr2,var_95_qr3,var_99_qr1,var_99_qr2,var_99_qr3)
write.csv(result_hs,"/Users/lisiqi/Desktop/result_hs.csv")

result_sh <- read.csv("/Users/lisiqi/Desktop/result_sh.csv")
result_hs <- read.csv("/Users/lisiqi/Desktop/result_hs.csv")
VaR95 <- abs(c(result_sh$var_95_qr1,result_sh$var_95_qr2,result_sh$var_95_qr3,result_hs$var_95_qr1,result_hs$var_95_qr2,result_hs$var_95_qr3))
t <- rep(c("SH.A1","SH.A2","SH.A3","HS1","HS2","HS3"),each=300)
data_VaR95 <- data.frame(VaR95,t)
boxplot(VaR95~t,data=data_VaR95,main="95% VaR",xlab="markets",ylab="VaR")

VaR99 <- abs(c(result_sh$var_99_qr1,result_sh$var_99_qr2,result_sh$var_99_qr3,result_hs$var_99_qr1,result_hs$var_99_qr2,result_hs$var_99_qr3))
data_VaR99 <- data.frame(VaR99,t)
boxplot(VaR99~t,data=data_VaR99,main="99% VaR",xlab="markets",ylab="VaR")