This section contains basic functions to support ‘HQGARCHEst’ and ‘HQGARCHRoll’.
#################################Load Packages###########################
library(quantreg) #qr
library(rugarch) #QMLE
###############################Functions#################################
h_update<-function(N,y,theta,intc) #iterative formula of h_t for GARCH(1,1)
{
#N is sample size; y is t.s. data from GARCH(1,1);
#theta is the parameter; (p+q+1)*1 vector
#intc is the initial value
alpha0=theta[1];alpha=theta[2];beta=theta[3]
sigma=NULL
sigma[1]=alpha0+alpha*intc+beta*intc #initial values for y_0^2 and h_0, both are intc
for (i in 2:N)
{
sigma[i]=alpha0+alpha*y[i-1]^2+beta*sigma[i-1]
}
return(sigma) #h_t
}
derivativeG<-function(N,x,theta,intc) #score vector and Hessian for Gaussian GARCH(1,1)
{
#N is sample size; orders p,q=1; x is N*1 vector
#theta is the QMLE parameter; (p+q+1)*1 vector
htilde=h_update(N,x,theta,intc) #N*1 vector
X1=c(intc,x[1:(N-1)]^2) #initial value for x_0^2 is intc
X2=c(intc,htilde[1:(N-1)]) #initial value for htilde_0 is intc
Zt=cbind(1,X1,X2) #N*(p+q+1) matrix
beta=theta[3]
P1h=matrix(0,N,3) #first derivative of h_t with respect to theta given initial values of P1h_0 being zero
P2h=array(0,dim=c(3,3,N)) #second derivative of h_t with respect to theta given initial values of P2h_0 being zero
P1ell=matrix(0,N,3) #first derivative of ell_t with respect to theta (score vector)
P2ell=array(0,dim=c(3,3,N)) #second derivative of ell_t with respect to theta (Hessian)
P1h[1,]=Zt[1,] #initial values of P1h_0 being zero
P1ell[1,]=(1-x[1]^2/htilde[1])/htilde[1]*P1h[1,]
P2ell[,,1]=(1-x[1]^2/htilde[1])/htilde[1]*P2h[,,1]+(2*x[1]^2/htilde[1]-1)/(htilde[1]^2)*tcrossprod(P1h[1,])
for (i in 2:N)
{
P1h[i,]=Zt[i,]+beta*P1h[i-1,]
P2h[,,i]=cbind(0,0,P1h[i-1,])+rbind(0,0,P1h[i-1,])+beta*P2h[,,i-1]
P1ell[i,]=(1-x[i]^2/htilde[i])/htilde[i]*P1h[i,]
P2ell[,,i]=(1-x[i]^2/htilde[i])/htilde[i]*P2h[,,i]+(2*x[i]^2/htilde[i]-1)/(htilde[i]^2)*tcrossprod(P1h[i,])
}
return(list(P1h=P1h,P2h=P2h,P1ell=P1ell,P2ell=P2ell))
}
QMLE_NR<-function(N,MM,p,q,tolerance,x,theta_int,intc) #QMLE for GARCH(1,1)
{
#N is sample size; orders p,q=1
#MM is maximum iteration
#tolerance is the convergence tolerance
#x is t.s. data; N*1 vector
#theta_int is initial values for parameter; (p+q+1)*1 vector
theta=matrix(0,p+q+1,MM)
i=1
laststep=1
theta[,i]<-theta_int
repeat
{
der=derivativeG(N,x,theta[,i],intc)
P1ells=apply(der$P1ell,2,sum) #summation of P1ell respect to time t from 1 to N
P2ells=apply(der$P2ell,c(1,2),sum) #summation of P2ell respect to time t from 1 to N
if (rcond(P2ells)<10^(-15))
{
#print("P2ells is computationally singular")
Indicator=1
break
}
theta[,i+1]=theta[,i]-solve(P2ells, tol = 1e-16)%*%P1ells
if (sum(theta[,i+1]>0)<3 | sum(theta[2:3,i+1])>2)
{
#print("The parameter has values being negative or greater than one")
Indicator=1
break
} #if the parameter has values being negative or greater than one, stop the repetition.
else
{
laststep=i+1
maxabs=max(abs(theta[,i+1]-theta[,i]))
i=i+1
#print(i)
if (i>=MM & maxabs>tolerance)
{
Indicator=2
break
} #if the iteration approaches the maximum iteration limitation but not converges, stop the repetition.
if (maxabs<=tolerance)
{
Indicator=0
break
} #if the iteration converges, stop the repetition.
}
}
return(list(Indicator=Indicator,theta_up=theta[,laststep],laststep=laststep))
}
psi<-function(x,tau)
{
tau-1*(x<0)
}
rk<-function(N,k,ehat0,ehat,tau,weight)
{
#tau is a scalar; k is the lag
#ehat0 are residuals from estimation; ehat are residuals or bootstrap residuals
#weight: (length(ehat)-k)*1 vector of equal weight in diagnosis or the random weight in bootstrap
var_a=var(abs(ehat0))
rk=sum(weight*psi(ehat[(k+1):N],tau)*abs(ehat[1:(N-k)]))/N/sqrt((tau-tau^2)*var_a)
return(rk)
}
CQEst<-function(tau,p=1,q=1,L,X,B)
{
#######################################################################
# This function is used to perform estimation based on the hybrid #
# estimation procedure and do diagnostic checking for conditional #
# quantiles based on the QACF of residuals #
#############################Input#####################################
# tau is the quantile level of interest #
# p and q are the orders in GARCH(p,q) model #
# L is a predetermined integer for the maximum lag in QACF #
# X is the t.s. data for analysis #
# B is size of Bootstrap #
#######################################################################
N=length(X) #N is the sample size
intc=mean(X[1:5]^2) #initial value of X_t^2 and h_t for t<=0
##############################Matrices of results#######################
thetahat<-matrix(0,p+q+1,1) #WCQE
R<-matrix(0,L,1)
QR<-NULL
thetatildeB<-matrix(0,p+q+1,1) #QMLE after correction in bootstrapping
thetahatB<-matrix(0,p+q+1,1) #WCQE in bootstrapping
dE<-matrix(0,p+q+1,B) #sqrt(N)*(thetahatB-thetahat)
Sigma1hatB<-matrix(0,p+q+1,p+q+1) #bootstrapped covariance matrices for WCQE
sdthetahatB<-matrix(0,p+q+1,1) #bootstrapped standard errors for WCQE
RB<-NULL #L*1 vector for bootstrap test statistics
dT<-matrix(0,L,B) #sqrt(N)*(RhatB-Rhat)
Sigma2hatB<-matrix(0,L,L) #bootstrapped covariance matrices for R
sdRB<-matrix(0,L,1) #bootstrapped standard errors for R
LCIRB<-UCIRB<-matrix(0,nrow=L,ncol=1) #confidence intervals for R
cvQ<-NULL #critical values for Q
##############################QMLE for GARCH(1,1)##########################
spec1_1=ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0),include.mean = FALSE),
distribution.model = "norm")
garch1_1=ugarchfit(spec = spec1_1, data = X, solver = 'hybrid', fit.control = list(stationarity = 1))
Theta_int=as.vector(coef(garch1_1)) #initial value for QMLE
QMLE=QMLE_NR(N,MM=20,p,q,tolerance=10^(-6),X,Theta_int,intc)
QMLEInd=QMLE$Indicator
if(QMLEInd!=0)
{
Ind=1#;print(Ind)
thetatilde=Theta_int #(p+q+1)*1 vector
}else
{
Ind=0#;print(Ind)
thetatilde=QMLE$theta_up #(p+q+1)*1 vector
}
######################GARCH(p,q) two-step estimation#######################
Y=X^2*sign(X) #Responose
#####step I. QMLE for GARCH(1,1) to obtain h_t
httilde=h_update(N,X,thetatilde,intc) #length(httilde)=N
#####step II. WCQE at specific tau
#Create Regressor
X1=c(intc,X[1:(N-1)]^2)
X2=c(intc,httilde[1:(N-1)])
XX=cbind(X1,X2) #Regressor without intercept in step II of GARCH(1,1) model; N*(p+q)
fit=rq(Y ~ XX, tau=tau, weights=1/httilde)
thetahat=as.vector(coef(fit))
names(thetahat)=c("alpha0tau","alpha1tau","beta1tau")
######################GARCH(p,q) diagnostic checking#######################
ehat=(Y-as.vector(cbind(1,XX)%*%thetahat))/httilde
for (k in 1:L)
{
R[k]=rk(N,k,ehat,ehat,tau,weight=rep(1,N-k))
}
#######################Correction matrices##################################
Der=derivativeG(N,X,thetatilde,intc)
P1h=Der$P1h #N*(p+q+1)
P1ell=Der$P1ell #N*(p+q+1)
Jtilde=crossprod(P1h/httilde)/N #(p+q+1)*(p+q+1)
Jinv=solve(Jtilde, tol = 1e-16)
##################################Bootstrap#################################
set.seed(12345)
W<-matrix(rexp(N*B, rate = 1),N,B) # standard exponential
b=1
while (b<=B)
{
##########################Weights generation###############################
weight=W[,b] #standard exponential
############################Correction terms###############################
e_c=Jinv%*%apply((weight-1)*P1ell,2,mean) #correction term
####step I. corrected QMLE for GARCH(1,1) to obtain h_tB for bootstrap
thetatildeB=thetatilde-e_c #corrected QMLE
httildeB=h_update(N,X,thetatildeB,intc)
####step II. WCQE at specific tau for bootstrap
X2B=c(intc,httildeB[1:(N-1)])
XXB=cbind(X1,X2B) #Regressor without intercept in step II for bootstrap
fitB=rq(Y ~ XXB, tau=tau, weights=weight/httilde)
thetahatB=as.vector(coef(fitB))
dE[,b]=sqrt(N)*(thetahatB-thetahat)
##############GARCH(p,q) diagnostic checking for bootstrap##############
ehatB=(Y-as.vector(cbind(1,XXB)%*%thetahatB))/httilde
for (k in 1:L)
{
RB[k]=rk(N,k,ehat,ehatB,tau,weight=weight[(k+1):N])
}
dT[,b]=sqrt(N)*(RB-R)
b=b+1
#print(b)
}
#QMLE
kappa=mean(X^4/httilde^2)
sdthetatilde=sqrt(diag((kappa-1)*Jinv/N))
#WCQE
sdthetahatB=apply(dE,1,sd)/sqrt(N)
#R
Sigma2hatB=var(t(dT))
sdRB=sqrt(diag(Sigma2hatB)/N) #apply(dT,1,sd)/sqrt(N)
#check test statistic individually
LCIRB=apply(dT,1,quantile,probs=0.025)
UCIRB=apply(dT,1,quantile,probs=0.975)
#check test statistic jointly
QR=N*t(R)%*%solve(Sigma2hatB)%*%R
return(list(thetatilde=thetatilde,sdQMLE=sdthetatilde,
thetahat=thetahat,sdE=sdthetahatB,
R=R,sdR=sdRB,QR=QR,
LCIRB=LCIRB,UCIRB=UCIRB)
)
#############################Output#######################################
# thetatilde and sdQMLE are estimates and standard errors of QMLE #
# thetahat and sdE are estimates and standard errors of hybrid estimator #
# R and sdR are estimates and standard errors of QACFs #
# QR is the value of portmanteau test statistic Q(L) at maximum lag L #
# LCIRB and UCIRB are lower and upper bounds for QACF #
##########################################################################
}
CICQE<-function(tau,p=1,q=1,B,a,b,M,data)
{
#######################################################################
# This package is used to forecast conditional quantiles for the data #
# and to construct their confidence intervals (CIs) #
##############################Input####################################
# tau is the quantile level of interest #
# p and q are the orders in GARCH(p,q) model #
# B is size of Bootstrap #
# a and b determine the starting and ending points respectively #
# M is the out-sample size of forecasting #
# data is the t.s. data to construct CIs for conditional quantiles #
# NN=b-M #start point of prediction in original data #
#######################################################################
#########################################################################
#matrices of results
Ind<-NULL #indicator for convergence of QMLE
thetaH<-thetaQ<-matrix(0,p+q+1,M) #hybrid and QMLE
QuantileH<-CIL<-CIU<-NULL #one-step ahead CQE and its CI
QuantileB<-NULL #one-step ahead CQE in bootstrap
NN=b-M #start point of prediction in original data
n0=max(p,q)
spec1_1=ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(p, q)),
mean.model = list(armaOrder = c(0, 0),include.mean = FALSE),
distribution.model = "norm")
for (l in 1:M)
{
##############################Data generation############################
Xlog=log(data$Y[a:(NN+l-1)])
X=diff(Xlog)
n=length(X) #sample size
##############################Estimation#################################
########################GARCH(1,1) hybrid estimation#####################
####step I. QMLE for GARCH(1,1) to obtain h_t
spec1_1=ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(p, q)),
mean.model = list(armaOrder = c(0, 0),include.mean = FALSE),
distribution.model = "norm")
garch1_1=ugarchfit(spec = spec1_1, data = X, solver = 'hybrid', fit.control = list(stationarity = 1))
Theta_int=as.vector(coef(garch1_1)) #initial value for QMLE
intc=mean(X[1:5]^2) #initial value of X_t^2 and h_t for t<=0
#intc=Theta_int[1]/(1-Theta_int[2]-Theta_int[3]) #initial value of X_t^2 for zero or negative t
QMLE=QMLE_NR(n,MM=20,p=1,q=1,tolerance=10^(-6),X,Theta_int,intc)
QMLEInd=QMLE$Indicator
if(QMLEInd!=0)
{
Ind[l]=1#;print(Ind)
thetatilde=Theta_int #(p+q+1)*1 vector
}else
{
Ind[l]=0#;print(Ind)
thetatilde=QMLE$theta_up #(p+q+1)*1 vector
}
httilde=h_update(n,X,thetatilde,intc) #length(httilde)=n
####step II. WCQE at specific tau
#Create Responose and Regressor
X1=c(intc,X[1:(n-1)]^2)
X2=c(intc,httilde[1:(n-1)])
XX=cbind(X1,X2) #Regressor without intercept in step II of GARCH(1,1) model; n*(p+q)
Y=X^2*sign(X) #Responose
fit=rq(Y ~ XX, tau=tau, weights=1/httilde)
thetaH[,l]=as.vector(coef(fit))
#############################One-step ahead CQE############################
QHY=as.numeric(crossprod(thetaH[,l],c(1,X[n]^2,httilde[n])))
QuantileH[l]=sqrt(abs(QHY))*sign(QHY)
################################Bootstraped CI#############################
#######################correction matrices##################################
Der=derivativeG(n,X,thetatilde,intc)
P1h=Der$P1h #n*(p+q+1)
P1ell=Der$P1ell #n*(p+q+1)
Jtilde=crossprod(P1h/httilde)/n #(p+q+1)*(p+q+1)
Jinv=solve(Jtilde, tol = 1e-16)
set.seed(12345)
W<-matrix(rexp(n*B, rate = 1),n,B) # standard exponential
b=1
while (b<=B)
{
##########################weights generation###############################
weight=W[,b] #standard exponential
############################Correction terms###############################
e_c=Jinv%*%apply((weight-1)*P1ell,2,mean) #correction term
####step I. corrected QMLE for GARCH(1,1) to obtain h_tB for bootstrap
thetatildeB=thetatilde-e_c #corrected QMLE
httildeB=h_update(n,X,thetatildeB,intc)
####step II. WCQE at specific tau for bootstrap
X2B=c(intc,httildeB[1:(n-1)])
XXB=cbind(X1,X2B) #Regressor without intercept in step II for bootstrap
fitB=rq(Y ~ XXB, tau=tau, weights=weight/httilde)
thetahatB=as.vector(coef(fitB))
QHYB=as.numeric(crossprod(thetahatB,c(1,X[n]^2,httildeB[n])))
QuantileB[b]=sqrt(abs(QHYB))*sign(QHYB)
b=b+1
#print(b)
}
CIL[l]=2*QuantileH[l]-quantile(QuantileB,probs=0.975)
CIU[l]=2*QuantileH[l]-quantile(QuantileB,probs=0.025)
print(l)
}
return(list(QuantileH=QuantileH,CIL=CIL,CIU=CIU))
#############################Output#######################################
# QuantileH is the rolling forecast of the conditional quantile #
# CIL and CIU are lower and upper bounds for rolling forecasts #
##########################################################################
}
#################################Ending#####################################`HQGARCHEst’ provides a function to perform the proposed hybrid conditional quantile estimation and diagnostic checking. You can use it to obtain model estimates, portmanteau test statistic, residual QACFs and residual QACF plot.
source("functions.R")
HQGARCH <- function(x, ...) UseMethod("HQGARCH")
HQGARCH.default <- function(tau=0.05,p=1,q=1,L=30,x,B=1000)
{
#############################Input#####################################
# tau is the quantile level of interest #
# p and q are the orders in GARCH(p,q) model #
# L is a predetermined integer for the maximum lag in QACF #
# x is the t.s. data for analysis #
# B is size of Bootstrap #
#######################################################################
tau <- as.numeric(tau)
p <- as.integer(p)
q <- as.integer(q)
L <- as.integer(L)
x <- as.numeric(x)
B <- as.integer(B)
est <- CQEst(tau,p=1,q=1,L,x,B)
thetatilde <- est$thetatilde
sdQMLE <- est$sdQMLE
N <- length(x)
intc <- mean(x[1:5]^2)
httilde <- h_update(N,x,thetatilde,intc)
thetahat <- est$thetahat
sdE <- est$sdE
QHY <- cbind(1,x^2,httilde)%*%thetahat
R <- est$R
sdR <- est$sdR
LCIRB <- est$LCIRB/sqrt(N)
UCIRB <- est$UCIRB/sqrt(N)
QR <- est$QR
test.pval <- pchisq(QR,df=L,lower.tail = FALSE)
QACFPlot <- plot(c(0,L), c(-0.10,0.10), type = "n", xlab = "Lag", ylab = "Residual QACF", ylim=c(-0.10,0.10),xaxt="n", yaxt="n")+
axis(1,at=seq(0,L,5))+
axis(2,at=seq(-0.10,0.10,0.05))+
abline(h=0)+
for (i in 1:L)
{
lines(c(i,i),c(0,R[i]))
}+
lines(x=c(1:L),y=LCIRB,type="l",lty="dotted",col="black")+
lines(x=c(1:L),y=UCIRB,type="l",lty="dotted",col="black")
est$QACFPlot <- QACFPlot
est$coefficients <- thetahat
est$fitted.values <- as.vector(sqrt(abs(QHY))*sign(QHY))
est$residuals <- x - est$fitted.values
est$test.pval <- test.pval
Q <- matrix(QR,ncol=1); colnames(Q) <- "Q.value"
pval <- matrix(test.pval,ncol=1); colnames(pval) <- "p.value"
est$test <- cbind(Q.value = Q, p.value = pval)
est$call <- match.call()
class(est) <- "HQGARCH"
est
}
print.HQGARCH <- function(x, ...)
{
#est <- CQEst(tau,p,q,L,x=x,B)
cat("Call:\n")
print(x$call)
cat("\nCoefficients:\n")
print(x$coefficients)
cat("\nPortmanteau test:\n")
printCoefmat(x$test, P.value=TRUE, has.Pvalue=TRUE)
}
summary.HQGARCH <- function(object, ...)
{
se <- object$sdE
tval <- coef(object) / se
EstTAB <- cbind(Estimate = coef(object),
StdErr = se,
t.value = tval,
p.value = 2*pnorm(abs(tval),lower.tail = FALSE))
QACF <- matrix(object$R,ncol=1); colnames(QACF) <- "QACF"
QACFTAB <- cbind(QACF = QACF,
StdErr = object$sdR,
Lower.CI = object$LCIRB,
Upper.CI = object$UCIRB)
QACFPlot <- object$QACFPlot
#Q <- matrix(object$QR,ncol=1); colnames(Q) <- "Q.value"
#pval <- matrix(object$test.pval,ncol=1); colnames(pval) <- "p.value"
TestTAB <- object$test #cbind(Q.value = Q, p.value = pval)
res <- list(call=object$call,
coefficients=EstTAB,
QACFs=QACFTAB,
QACFPlot=QACFPlot,
test=TestTAB)
class(res) <- "summary.HQGARCH"
res
}
print.summary.HQGARCH <- function(x, ...)
{
cat("Call:\n")
print(x$call)
cat("\nCoefficients:\n")
printCoefmat(x$coefficients, P.value=TRUE, has.Pvalue=TRUE)
cat("\nPortmanteau test:\n")
printCoefmat(x$test, P.value=TRUE, has.Pvalue=TRUE)
cat("\nResidual QACFs:\n")
printCoefmat(x$QACFs)
cat("\nResidual QACF Plot:\n")
print(x$QACFPlot)
}
#################################Ending####################################‘HQGARCHRoll’ is used to do rolling forecasting for conditional quantiles. We can obtain rolling forecasts and their 95% confidence intervals.
source("functions.R")
HQGARCHRoll <- function(data, ...) UseMethod("HQGARCHRoll")
HQGARCHRoll.default <- function(tau=0.05,p=1,q=1,B=1000,a,b,M,data)
{
##############################Input####################################
# tau is the quantile level of interest #
# p and q are the orders in GARCH(p,q) model #
# B is size of Bootstrap #
# a and b determine the starting and ending points respectively #
# M is the out-sample size of forecasting #
# data is the t.s. data to construct CIs for conditional quantiles #
#######################################################################
tau <- as.numeric(tau)
p <- as.integer(p)
q <- as.integer(q)
B <- as.integer(B)
a <- as.integer(a)
b <- as.integer(b)
M <- as.integer(M)
data <- as.data.frame(data)
rollingfc <- CICQE(tau,p=1,q=1,B,a,b,M,data)
QuantileH <- rollingfc$QuantileH
CIL <- rollingfc$CIL
CIU <- rollingfc$CIU
X <- diff(log(data$Y[a:b]))
N <- length(X)
RollPlot <- plot(X,type='l',ylab="",xlab="",cex.axis=1,ylim=c(-0.15, 0.15))+
lines(x=c((N-M+1):N),y=QuantileH,type="l",lty="solid",lwd=0.1,col="blue")+
lines(x=c((N-M+1):N),y=CIL,type="l",lty="solid",lwd=0.1,col="red")+
lines(x=c((N-M+1):N),y=CIU,type="l",lty="solid",lwd=0.1,col="red")+
mtext("Log return", side = 2, line = 2.5, cex=1)+
mtext("Time", side = 1, line = 2.5, cex=1)
rollingfc$RollPlot <- RollPlot
rollingfc$call <- match.call()
class(rollingfc) <- "HQGARCHRoll"
rollingfc
}
print.HQGARCHRoll <- function(x, ...)
{
#est <- CQEst(tau,p,q,L,x=x,B)
cat("Call:\n")
print(x$call)
cat("\nRolling forecasting:\n")
print(x$RollPlot)
}
#################################Ending####################################This section provides R codes for the empirical data analysis of the paper, which illustrates how to use ‘HQGARCHEst’ and ‘HQGARCHRoll’.
setwd("C:/Users/admin/Desktop")
source("functions.R")
source("HQGARCHEst.R")
source("HQGARCHRoll.R")
#################################S&P500#####################################
#load data
data <- read.table("SP500.txt", header=TRUE) #consider S&P500 as an example
attach(data)
Time <- strptime(data$Date, format = "%m/%d/%Y")
SPI <- as.numeric(data$Y)
start <- strptime("2008-01-02",format = "%Y-%m-%d") #starting point
end <- strptime("2016-06-30",format = "%Y-%m-%d") #ending point
Fstart <- strptime("2010-01-04",format = "%Y-%m-%d")
Ind <- (1:length(SPI))[Time>=start & Time<=end]
t <- Time[Ind]
Y0 <- SPI[Ind]
X <- diff(log(Y0)) #log return
N <- length(X) #sample size
a <- (1:length(SPI))[Time==start] #data origin
b <- (1:length(SPI))[Time==end] #data ending point
M <- length(Ind[Time>=Fstart & Time<=end]) #out-sample size of forecasting
################################Settings###################################
tau <- 0.05 #specific quantile
B <- 1000 #size of Bootstrap
p <- 1;q <- 1 #order in GARCH(p,q)
L <- 30 #maximum lag in QACF
################################Analysis###################################
#estimation and diagnotic checking
EstResult <- HQGARCH(tau=0.05,p=1,q=1,L=30,x=X,B=1000)
EstResult
summary(EstResult)
#rolling forecasting
RollResult <- HQGARCHRoll(tau=0.05,p=1,q=1,B=1000,a=1759,b=3898,M=1635,data=data)
RollResult
#################################Ending####################################