Nonparametric volatility prediction

Introduction

This is the home page of the article "Nonparametric volatility prediction", written by Jussi Klemelä.

Abstract of the article

Volatility prediction can be defined as the estimation of the conditional expectation of a squared return of a financial asset. Many of the classical nonparametric regression estimators can be applied in volatility prediction. Examples of nonparametric estimators include moving averages and kernel estimators. However, it has been difficult to beat some parametric estimators from the GARCH family using nonparametric estimators. We review some promising suggestions for nonparamteric volatility prediction.

Reference

@article{wiresvola,
title = {Nonparametric volatility prediction},
author = {Klemel\"a, J.},
year = {2020},
journal = {{WIRE}s Comput. Statist.},
volume = {12},
number = {3},
note = {e1491},
web = {https://doi.org/10.1002/wics.1491},
}

Reproducing the computations

Install software

source("http://jklm.fi/regpro/regpro.R")

Read SP500 Data

file<-"http://jklm.fi/statopti/GSPC.csv"
gspc0<-read.csv(file=file)
gspc<-gspc0[,6]

n0<-length(gspc)
retu<-gspc[2:n0]/gspc[1:(n0-1)]    
r<-log(retu)
n<-length(r)

times<-matrix(0,n,1)
delta<-1/251.63
alku<-1950+0/12
for (i in 1:length(times)) times[i]<-alku+(i-1)*delta

Figure 1

t0<-1000

# GARCH(1,1) for SP500: Rolling

library("tseries")

roll<-2500 
t0<-250
step<-1
ota<-seq(t0+1,n,step)

alpha0t<-matrix(0,n,1)
alpha1t<-matrix(0,n,1)
betat<-matrix(0,n,1)
sigmat<-matrix(0,n,1)
for (t in 1:t0){
   now<-r[1:t]
   sigmat[t]<-sd(now)
}
for (i in 1:length(ota)){
    t<-ota[i]
    #
    if (i==length(ota)) up<-n else up<-ota[i+1]-1
    if (is.null(roll)) beg<-1 else beg<-max(1,t-roll)
    now<-r[beg:t]
    #
    ga<-garch(now,control=garch.control(trace=FALSE))   #-mean(y))  
    alpha0<-as.numeric(ga$coef[1])  
    alpha1<-as.numeric(ga$coef[2])  
    beta<-as.numeric(ga$coef[3])    
    alpha0t[t]<-alpha0
    alpha1t[t]<-alpha1
    betat[t]<-beta
    sigmat[t]<-sqrt(alpha0+alpha1*r[t]^2+beta*sigmat[t-1]^2)
}
sigmat[1]<-sigmat[2]

# one-day prediction

hori<-1
     
potti<-matrix(0,n,2)
potti[,1]<-r^2
for (i in (hori+1):n) potti[i,2]<-sigmat[i-hori]^2

# prediction of realized volatility

realvola<-matrix(0,n,1)
observed<-matrix(0,n,1)
hori<-20
for (i in (1+hori):n){
    observed[i]<-sum(r[(i-hori+1):i]^2)
    #
    for (eta in 1:hori){
	 sigma<-sigmat[i-hori]
         a0<-alpha0t[i-hori] 
         a1<-alpha1t[i-hori] 
         b<-betat[i-hori]    
         barsigma2<-a0/(1-a1-b)
         vest<-barsigma2+(a1+b)^(eta-1)*(sigma^2-barsigma2)
         #
         realvola[i]<-realvola[i]+vest	       
    }	   	       
} 
sotti<-matrix(0,n,2)
sotti[,1]<-observed
sotti[,2]<-realvola

# Figure 1

cex.axis<-1.5
cex.lab<-1.5
cex.sub<-1.5
cex<-1.5
lty<-1.5

ota<-c(1000:n)
matplot(times[ota],potti[ota,],type="l",lty=1,ylim=c(0,0.005),
xlab="",ylab="",cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
mtext("squared return",side=2,cex=1.5,line=2.4)
title(sub="(a)",cex.sub=cex.sub)

ota<-c((n-400):n)
matplot(times[ota],potti[ota,],type="l",lty=1,ylim=c(0,0.0006),
xlab="",ylab="",cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
mtext("squared return",side=2,cex=1.5,line=2.4)
title(sub="(b)",cex.sub=cex.sub)

ota<-c(1000:n)
matplot(times[ota],sotti[ota,],type="l",lty=1,ylim=c(0,0.008),
xlab="",ylab="",cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
mtext("realized volatility",side=2,cex=1.5,line=2.4)
title(sub="(c)",cex.sub=cex.sub)

ota<-c((n-400):n)
matplot(times[ota],sotti[ota,],type="l",lty=1, #,ylim=c(0,0.008))
xlab="",ylab="",cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
mtext("realized volatility",side=2,cex=1.5,line=2.4)
title(sub="(d)",cex.sub=cex.sub)

Figure 2

t0<-1000
eta<-1

# GARCH(1,1) for SP500: Rolling

library("tseries")

roll<-2500 
t0<-250
step<-1
ota<-seq(t0+1,n,step)

alpha0t<-matrix(0,n,1)
alpha1t<-matrix(0,n,1)
betat<-matrix(0,n,1)
sigmat<-matrix(0,n,1)
for (t in 1:t0){
   now<-r[1:t]
   sigmat[t]<-sd(now)
}
for (i in 1:length(ota)){
    t<-ota[i]
    #
    if (i==length(ota)) up<-n else up<-ota[i+1]-1
    if (is.null(roll)) beg<-1 else beg<-max(1,t-roll)
    now<-r[beg:t]
    #
    ga<-garch(now,control=garch.control(trace=FALSE))   #-mean(y))  
    alpha0<-as.numeric(ga$coef[1])  
    alpha1<-as.numeric(ga$coef[2])  
    beta<-as.numeric(ga$coef[3])    
    alpha0t[t]<-alpha0
    alpha1t[t]<-alpha1
    betat[t]<-beta
    sigmat[t]<-sqrt(alpha0+alpha1*r[t]^2+beta*sigmat[t-1]^2)
}
sigmat[1]<-sigmat[2]
     
gcspe<-matrix(0,n,1)
diffet<-matrix(0,n,1)
for (i in (t0+1):(n-eta)){
    sigma<-sigmat[i]
    a0<-alpha0t[i] 
    a1<-alpha1t[i] 
    b<-betat[i]    
    barsigma2<-a0/(1-a1-b)
    hava<-r[i+eta]
    vest<-barsigma2+(a1+b)^(eta-1)*(sigma^2-barsigma2)
    diffet[i+eta]<-abs(vest-hava^2)^2
    gcspe[i+eta]<-sum(diffet[1:(i+eta)],na.rm=TRUE)
}

# Heston-Nandi for SP500: Rolling

library(fOptions)

roll<-2500
t0<-250
alpha0t<-matrix(0,n,1)
alpha1t<-matrix(0,n,1)
betat<-matrix(0,n,1)
gammat<-matrix(0,n,1)
sigmat<-matrix(0,n,1)
step<-20
ota<-seq(t0+1,n,step)
for (t in 1:t0){
   now<-r[1:t]
   sigmat[t]<-sd(now)
}
for (i in 1:length(ota)){
    t<-ota[i]
    if (i==length(ota)) up<-n else up<-ota[i+1]-1
    if (is.null(roll)) beg<-1 else beg<-max(1,t-roll)
    now<-r[beg:t]
    #
    hn<-hngarchFit(now)
    alpha0<-hn$model$omega
    alpha1<-hn$model$alpha  
    beta<-hn$model$beta
    gamma<-hn$model$lambda  
    alpha0t[t:up]<-alpha0
    alpha1t[t:up]<-alpha1
    betat[t:up]<-beta
    gammat[t:up]<-gamma
}
for (t in (t0+1):length(sigmat)){
    alpha0<-alpha0t[t]
    alpha1<-alpha1t[t]  
    beta<-betat[t]
    gamma<-gammat[t]  
    sigmat[t]<-
sqrt(alpha0+alpha1*(r[t]/sigmat[t-1]-gamma*sigmat[t-1])^2+beta*sigmat[t-1]^2)
}
sigmat[1]<-sigmat[2]

hn.cspe<-matrix(0,n,1)
diffet<-matrix(0,n,1)
for (i in (t0+1):(n-eta)){
    sigma<-sigmat[i]
    hava<-r[i+eta]
    vest<-sigma^2
    diffet[i+eta]<-abs(vest-hava^2)^2
    hn.cspe[i+eta]<-sum(diffet[1:(i+eta)],na.rm=TRUE)
}

# asymmetric MA cross-validation #######
lambdat<-c(0,0.1,0.2,0,3)
hoot<-c(5,10,20,40)
gammat<-exp(-1/hoot)
acspet<-matrix(0,n,length(hoot)*length(lambdat))
aspet<-matrix(0,n,length(hoot)*length(lambdat))
for (ll in 1:length(lambdat)){
for (gg in 1:length(gammat)){
      lambda<-lambdat[ll]
      gamma<-gammat[gg]
      diffet<-matrix(0,n,1)
      val<-mean(r[1:t0]^2)
      for (t in (t0+1):(n-eta)){
         hava<-r[t+eta]
         val<-(1-gamma)*val*(r[t]/sqrt(val)-lambda)^2+gamma*val
         diffet[t+eta]<-abs(val-hava^2)^2
         #
         ind<-length(gammat)*(ll-1)+gg 			
         acspet[t+eta,ind]<-sum(diffet[1:(t+eta)],na.rm=TRUE)
         aspet[t+eta,ind]<-abs(val-hava^2)^2
      }
}}
acspe.cv<-matrix(0,n,1)
for (t in (t0+1):(n-eta)){
    ine<-which.min(acspet[t,])
    acspe.cv[t]<-acspe.cv[t-1]+aspet[t,ine]
}

# Figure 2
		 
cex.axis<-1.5
cex.lab<-1.5
cex.sub<-1.5
cex<-1.5
lty<-1.5

ero1<-hn.cspe-gcspe     

ero2<-acspe.cv-gcspe           
		 
ota<-c(t0:(n-100))
plot(times[ota],ero1[ota],type="l",lty=1,#ylim=c(-0.0001,0.0001),
xlab="",ylab="",cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
#text(2000,0.00014,"Heston-Nandi",col="black",cex=cex)
segments(1,0,30000,0,col="green")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)
title(sub="(a)",cex.sub=cex.sub)

ota<-c(t0:(n-100))
plot(times[ota],ero2[ota],type="l",lty=1,#ylim=c(-0.0001,0.0001),
xlab="",ylab="",cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
#text(2000,0.00014,"EWMA",col="black",cex=cex)
segments(1,0,30000,0,col="green")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)
title(sub="(b)",cex.sub=cex.sub)

June 2019