Introduction

This is the home page of the article "Volatility Prediction Using Kernel Regression", written by Jussi Klemelä.

The article is available in volapred.pdf.

Supplementary material is available in volapred-supple.pdf.

Supplementary material about Dow Jones 30 stocks is available in volapred-supple-dj30.pdf.

Abstract of the article

We study the prediction of a squared return of a financial asset in a discrete time setting. The predictor is a kernel regression estimator whose explanatory variables are a moving average of squared returns and a moving average of returns. The predictor performs better than the GARCH($1,1$) predictor and some related predictors, in the sense of the mean squared prediction error, and when applied in the estimation of conditional quantiles. The kernel predictor is able to cope with the leverage effect.

Reproducing the computations

Install software

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

Read Daily SP500 Data

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
times0<-matrix(0,length(y),1)
delta<-1/251.63
alku<-1950+2/12
for (i in 1:length(times0)) times0[i]<-alku+(i-1)*delta
plot(times0,y,type="l")
  
m<-length(y)
dendat<-y[2:m]/y[1:(m-1)]
times<-times0[2:m]
plot(times,dendat,type="l")

Scatter plots

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
times0<-matrix(0,length(y),1)
delta<-1/251.63
alku<-1950+2/12
for (i in 1:length(times0)) times0[i]<-alku+(i-1)*delta

m<-length(y)	
r<-y[2:m]/y[1:(m-1)]-1
times<-times0[2:m]
n<-length(r)

h2<-40
x00<-matrix(0,n,2)
for (i in 1:n){ 
   now<-r[1:i]
   x00[i,1]<-sqrt(ma(now^2,h=h2))
   x00[i,2]<-ma(now,h=h2)
}
x0<-copula.trans(x00)

plot(x00,xlab="",ylab="")
mtext("MA of squared return",side=1,cex=1.5,line=2.4)
mtext("MA of return",side=2,cex=1.5,line=2.4)

plot(x0,xlab="",ylab="")
mtext("MA of squared return",side=1,cex=1.5,line=2.4)
mtext("MA of return",side=2,cex=1.5,line=2.4)

GARCH estimation

m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)

library("tseries")   
t0<-250
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 (t in (t0+1):length(sigmat)){
    now<-r[1: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)
}

Asymmetric GARCH estimation (Heston-Nandi model)

m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)

library(fOptions)

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
    now<-r[1: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(asigmat)){
    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)
}

Read GARCH(1,1) Volatility

file<-"http://jussiklemela.com/statfina/sigmat"
sigmat<-as.matrix(read.csv(file=file))

times<-matrix(0,length(sigmat),1)
delta<-1/251.63
alku<-1950+2/12
for (i in 1:length(times)) times[i]<-alku+(i-1)*delta
plot(times,sigmat,type="l")

Read Asymmetric GARCH(1,1) Volatility

file<-"http://jussiklemela.com/statfina/sigmat-hn"
sigmat<-as.matrix(read.csv(file=file))

times<-matrix(0,length(sigmat),1)
delta<-1/251.63
alku<-1950+2/12
for (i in 1:length(times)) times[i]<-alku+(i-1)*delta
plot(times,sigmat,type="l")

Plot GARCH(1,1) Parameter Estimates

file<-"http://jussiklemela.com/statfina/garch-netret.csv"
para<-as.matrix(read.csv(file=file))
betat<-para[,4]
       
file<-"http://jussiklemela.com/statfina/garch-netret-hn.csv"
para<-as.matrix(read.csv(file=file))
h.betat<-para[,4]
gammat<-para[,5]
       
times<-matrix(0,length(betat),1)
delta<-1/251.63
alku<-1950+2/12
for (i in 1:length(times)) times[i]<-alku+(i-1)*delta

ota<-c(2000:(n-10))
plot(times[ota],betat[ota],type="l",xlab="",ylab="",
ylim=c(min(betat[ota],h.betat[ota]),max(betat[ota],h.betat[ota])))
matplot(times[ota],h.betat[ota],col="red",add=TRUE,type="l")
mtext(expression(beta),side=2,cex=1.5,line=2.4)

ota<-c(2000:(n-10))
plot(times[ota],gammat[ota],type="l",xlab="",ylab="",col="red")
mtext(expression(lambda),side=2,cex=1.5,line=2.4)

Comparison between (asymmetric) GARCH and (asymmetric) moving average

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)

t0<-250    

##########################
eta<-1

# GARCH
file<-"http://jussiklemela.com/statfina/sigmat"
sigmat<-as.matrix(read.csv(file=file))
#
gcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,1)
for (i in (t0+1):(n-eta)){
     sigma<-sigmat[i]
     hava<-r[i+eta]
     vest<-sigma^2
     diffet[i]<-abs(vest-hava^2)^2
     gcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}

# asymmetric GARCH
file<-"http://jussiklemela.com/statfina/sigmat-hn"
sigmat<-as.matrix(read.csv(file=file))
#
hcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,1)
for (i in (t0+1):(n-eta)){
     sigma<-sigmat[i]
     hava<-r[i+eta]
     vest<-sigma^2
     diffet[i]<-abs(vest-hava^2)^2
     hcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}

# asymmetric MA 
t0<-250
lambdat<-c(0,0.2)
h<-20
gamma<-exp(-1/h)
cspe1<-matrix(0,n,length(lambdat))
for (ll in 1:length(lambdat)){
      lambda<-lambdat[ll]
      diffet<-matrix(0,n-eta,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]<-abs(val-hava^2)^2
         cspe1[t,ll]<-sum(diffet[1:t],na.rm=TRUE)
      }
}

# asymmetric MA with cross-validation
lambdat<-c(0,0.2)
hoot<-c(5,10,20,30,40)
cspem1<-matrix(0,n,length(hoot))
cspem2<-matrix(0,n,length(hoot))
spem1<-matrix(0,n,length(hoot))
spem2<-matrix(0,n,length(hoot))
for (hh in 1:length(hoot)){
      h<-hoot[hh]
      gamma<-exp(-1/h)
      #
      lambda<-lambdat[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
         spem1[t,hh]<-abs(val-hava^2)^2
         cspem1[t,hh]<-sum(spem1[1:t,hh],na.rm=TRUE)
      }
      lambda<-lambdat[2]
      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
         spem2[t,hh]<-abs(val-hava^2)^2
         cspem2[t,hh]<-sum(spem2[1:t,hh],na.rm=TRUE)
      }
}
#
cspem11<-matrix(0,n,1)
for (i in t0:n){
    ine<-which.min(cspem1[i-eta,])
    cspem11[i]<-cspem11[i-1]+spem1[i,ine]
}
cspem22<-matrix(0,n,1)
for (i in t0:n){
    ine<-which.min(cspem2[i-eta,])
    cspem22[i]<-cspem22[i-1]+spem2[i,ine]
}

# asymmetric MA with aggregation
lambdat<-c(0,0.2)
hoot<-c(5,10,20,30,40)
cspem1<-matrix(0,n,length(hoot))
cspem2<-matrix(0,n,length(hoot))
spem1<-matrix(0,n,length(hoot))
spem2<-matrix(0,n,length(hoot))
pred1<-matrix(0,n,length(hoot))
pred2<-matrix(0,n,length(hoot))
for (hh in 1:length(hoot)){
      h<-hoot[hh]
      gamma<-exp(-1/h)
      #
      lambda<-lambdat[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
         spem1[t,hh]<-abs(val-hava^2)^2
         cspem1[t,hh]<-sum(spem1[1:t,hh],na.rm=TRUE)
	 pred1[t,hh]<-val   
      }
      lambda<-lambdat[2]
      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
         spem2[t,hh]<-abs(val-hava^2)^2
         cspem2[t,hh]<-sum(spem2[1:t,hh],na.rm=TRUE)
	 pred2[t,hh]<-val   		
      }
}
#
cspem11<-matrix(0,n,1)
minhootm11<-matrix(0,n,1)
aggre1<-matrix(0,n,1)
painot1<-matrix(0,n,length(hoot))	       
for (i in (t0+2):n){
    ine<-which.min(cspem1[i-1,])
    cspem11[i]<-cspem11[i-1]+spem1[i,ine]
    minhootm11[i]<-hoot[ine]
    #
    paino<-exp(-cspem1[i-1,])/sum(exp(-cspem1[i-1,]))
    aggre1[i]<-sum(pred1[i,]*paino)
    painot1[i,]<-paino
}
cspe1a<-matrix(0,n,1)
spe1a<-matrix(0,n,1)
for (t in (t0+1):(n-eta)){
    hava<-r[t+eta]
    val<-aggre1[t]
    spe1a[t]<-abs(val-hava^2)^2
    cspe1a[t]<-sum(spe1a[1:t],na.rm=TRUE)      
}
cspem22<-matrix(0,n,1)
minhootm22<-matrix(0,n,1)
aggre2<-matrix(0,n,1)
painot2<-matrix(0,n,length(hoot))
for (i in t0:n){
    ine<-which.min(cspem2[i-1,])
    cspem22[i]<-cspem22[i-1]+spem2[i,ine]
    minhootm22[i]<-hoot[ine]
    #
    paino<-exp(-cspem2[i-1,])/sum(exp(-cspem2[i-1,]))
    aggre2[i]<-sum(pred2[i,]*paino)
    painot2[i,]<-paino
}
cspe2a<-matrix(0,n,1)
spe2a<-matrix(0,n,1)
for (t in (t0+1):(n-eta)){
    hava<-r[t+eta]
    val<-aggre2[t]
    spe2a[t]<-abs(val-hava^2)^2
    cspe2a[t]<-sum(spe2a[1:t],na.rm=TRUE)      
}


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

ero<-hcspe-gcspe
ero11<-cspem11-gcspe
ero22<-cspem22-gcspe
ero33<-cspe1a-gcspe
ero44<-cspe2a-gcspe
#      
ota<-c(1000:(n-100))
plot(times[ota],ero[ota],type="l",xlab="",ylab="",
ylim=c(min(ero[ota],ero11[ota],ero22[ota],ero33[ota],ero44[ota]),
max(ero[ota],ero11[ota],ero22[ota],ero33[ota],ero44[ota])),
lty=1)
matplot(times[ota],ero11[ota],col="red",add=TRUE,type="l",lty=1)
matplot(times[ota],ero22[ota],col="blue",add=TRUE,type="l",lty=1)
matplot(times[ota],ero33[ota],col="orange",add=TRUE,type="l",lty=1)
matplot(times[ota],ero44[ota],col="violet",add=TRUE,type="l",lty=1)
segments(1,0,30000,0,col="green")

ota<-c(1000:(n-6500))
plot(times[ota],ero[ota],type="l",xlab="",ylab="",lty=1,
ylim=c(-0.4e-05,0.5e-05))
colos<-c("red","blue","orange","violet")
matplot(times[ota],ero11[ota,],col=colos[1],add=TRUE,type="l",lty=1)
matplot(times[ota],ero22[ota,],col=colos[2],add=TRUE,type="l",lty=1)
matplot(times[ota],ero33[ota,],col=colos[3],add=TRUE,type="l",lty=1)
matplot(times[ota],ero44[ota,],col=colos[4],add=TRUE,type="l",lty=1)
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)
segments(1,0,30000,0,col="green")

ota<-c(9000:(n-100))
matplot(times[ota],ero[ota],type="l",xlab="",ylab="",lty=1,
ylim=c(-4.0e-05,12.0e-05))
colos<-c("red","blue","orange","violet")
matplot(times[ota],ero11[ota,],col=colos[1],add=TRUE,type="l",lty=1)
matplot(times[ota],ero22[ota,],col=colos[2],add=TRUE,type="l",lty=1)
matplot(times[ota],ero33[ota,],col=colos[3],add=TRUE,type="l",lty=1)
matplot(times[ota],ero44[ota,],col=colos[4],add=TRUE,type="l",lty=1)
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)
segments(1,0,30000,0,col="green")
      
##########################
eta<-10

# GARCH
file<-"http://jussiklemela.com/statfina/garch-netret.csv"
para<-as.matrix(read.csv(file=file))
sigmat<-para[,1]
alpha0t<-para[,2]
alpha1t<-para[,3]		   
betat<-para[,4]
#
gcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,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)
    vest<-barsigma2+(a1+b)^(eta-1)*(sigma^2-barsigma2)
    hava<-r[i+eta]
    diffet[i]<-abs(vest-hava^2)^2
    gcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}

# asymmetric GARCH
file<-"http://jussiklemela.com/statfina/garch-netret-hn.csv"
para<-as.matrix(read.csv(file=file))
sigmat<-para[,1]
alpha0t<-para[,2]
alpha1t<-para[,3]		   
betat<-para[,4]
gammat<-para[,5]
#
hcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,1)
for (i in (t0+1):(n-eta)){
    sigma<-sigmat[i]
    a0<-alpha0t[i] 
    a1<-alpha1t[i]
    b<-betat[i]
    gamma<-gammat[i]
    barsigma2<-(a0+a1)/(1-a1*gamma^2-b)
    vest<-barsigma2+(a1*gamma^2+b)^(eta-1)*(sigma^2-barsigma2)   
    hava<-r[i+eta]
    diffet[i]<-abs(vest-hava^2)^2
    hcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}

# asymmetric MA 
t0<-250
lambdat<-c(0,0.2)
h<-40
gamma<-exp(-1/h)
cspe1<-matrix(0,n,length(lambdat))
for (ll in 1:length(lambdat)){
      lambda<-lambdat[ll]
      diffet<-matrix(0,n-eta,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]<-abs(val-hava^2)^2
         cspe1[t,ll]<-sum(diffet[1:t],na.rm=TRUE)
      }
}

# asymmetric MA with cross-validation
lambdat<-c(0,0.2)
hoot<-c(20,30,40,60,80,100,120,140,200,300)
cspem1<-matrix(0,n,length(hoot))
cspem2<-matrix(0,n,length(hoot))
spem1<-matrix(0,n,length(hoot))
spem2<-matrix(0,n,length(hoot))
for (hh in 1:length(hoot)){
      h<-hoot[hh]
      gamma<-exp(-1/h)
      #
      lambda<-lambdat[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
         spem1[t,hh]<-abs(val-hava^2)^2
         cspem1[t,hh]<-sum(spem1[1:t,hh],na.rm=TRUE)
      }
      lambda<-lambdat[2]
      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
         spem2[t,hh]<-abs(val-hava^2)^2
         cspem2[t,hh]<-sum(spem2[1:t,hh],na.rm=TRUE)
      }
}
#
cspem11<-matrix(0,n,1)
for (i in t0:n){
    ine<-which.min(cspem1[i-eta,])
    cspem11[i]<-cspem11[i-1]+spem1[i,ine]
}
cspem22<-matrix(0,n,1)
for (i in t0:n){
    ine<-which.min(cspem2[i-eta,])
    cspem22[i]<-cspem22[i-1]+spem2[i,ine]
}

# asymmetric MA with aggregation
lambdat<-c(0,0.2)
hoot<-c(20,30,40,60,80,100,120,140,200,300)
cspem1<-matrix(0,n,length(hoot))
cspem2<-matrix(0,n,length(hoot))
spem1<-matrix(0,n,length(hoot))
spem2<-matrix(0,n,length(hoot))
pred1<-matrix(0,n,length(hoot))
pred2<-matrix(0,n,length(hoot))
for (hh in 1:length(hoot)){
      h<-hoot[hh]
      gamma<-exp(-1/h)
      #
      lambda<-lambdat[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
         spem1[t,hh]<-abs(val-hava^2)^2
         cspem1[t,hh]<-sum(spem1[1:t,hh],na.rm=TRUE)
	 pred1[t,hh]<-val   
      }
      lambda<-lambdat[2]
      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
         spem2[t,hh]<-abs(val-hava^2)^2
         cspem2[t,hh]<-sum(spem2[1:t,hh],na.rm=TRUE)
	 pred2[t,hh]<-val   		
      }
}
#
cspem11<-matrix(0,n,1)
minhootm11<-matrix(0,n,1)
aggre1<-matrix(0,n,1)
painot1<-matrix(0,n,length(hoot))	       
for (i in (t0+2):n){
    ine<-which.min(cspem1[i-eta,])
    cspem11[i]<-cspem11[i-1]+spem1[i,ine]
    minhootm11[i]<-hoot[ine]
    #
    paino<-exp(-cspem1[i-eta,])/sum(exp(-cspem1[i-eta,]))
    aggre1[i]<-sum(pred1[i,]*paino)
    painot1[i,]<-paino
}
cspe1a<-matrix(0,n,1)
spe1a<-matrix(0,n,1)
for (t in (t0+1):(n-eta)){
    hava<-r[t+eta]
    val<-aggre1[t]
    spe1a[t]<-abs(val-hava^2)^2
    cspe1a[t]<-sum(spe1a[1:t],na.rm=TRUE)      
}
cspem22<-matrix(0,n,1)
minhootm22<-matrix(0,n,1)
aggre2<-matrix(0,n,1)
painot2<-matrix(0,n,length(hoot))
for (i in t0:n){
    ine<-which.min(cspem2[i-eta,])
    cspem22[i]<-cspem22[i-1]+spem2[i,ine]
    minhootm22[i]<-hoot[ine]
    #
    paino<-exp(-cspem2[i-eta,])/sum(exp(-cspem2[i-eta,]))
    aggre2[i]<-sum(pred2[i,]*paino)
    painot2[i,]<-paino
}
cspe2a<-matrix(0,n,1)
spe2a<-matrix(0,n,1)
for (t in (t0+1):(n-eta)){
    hava<-r[t+eta]
    val<-aggre2[t]
    spe2a[t]<-abs(val-hava^2)^2
    cspe2a[t]<-sum(spe2a[1:t],na.rm=TRUE)      
}


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

ero<-hcspe-gcspe
ero11<-cspem11-gcspe
ero22<-cspem22-gcspe
ero33<-cspe1a-gcspe
ero44<-cspe2a-gcspe
#      
ota<-c(1000:(n-100))
plot(times[ota],ero[ota],type="l",xlab="",ylab="",
ylim=c(min(ero[ota],ero11[ota],ero22[ota],ero33[ota],ero44[ota]),
max(ero[ota],ero11[ota],ero22[ota],ero33[ota],ero44[ota])),
lty=1)
matplot(times[ota],ero11[ota],col="red",add=TRUE,type="l",lty=1)
matplot(times[ota],ero22[ota],col="blue",add=TRUE,type="l",lty=1)
matplot(times[ota],ero33[ota],col="orange",add=TRUE,type="l",lty=1)
matplot(times[ota],ero44[ota],col="violet",add=TRUE,type="l",lty=1)
segments(1,0,30000,0,col="green")

ota<-c(1000:(n-6500))
plot(times[ota],ero[ota],type="l",xlab="",ylab="",lty=1,
ylim=c(-0.2e-05,0.9e-05))
colos<-c("red","blue","orange","violet")
matplot(times[ota],ero11[ota,],col=colos[1],add=TRUE,type="l",lty=1)
matplot(times[ota],ero22[ota,],col=colos[2],add=TRUE,type="l",lty=1)
matplot(times[ota],ero33[ota,],col=colos[3],add=TRUE,type="l",lty=1)
matplot(times[ota],ero44[ota,],col=colos[4],add=TRUE,type="l",lty=1)
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)
segments(1,0,30000,0,col="green")

ota<-c(9000:(n-100))
matplot(times[ota],ero[ota],type="l",xlab="",ylab="",lty=1,
ylim=c(-14.0e-05,4.0e-05))
colos<-c("red","blue","orange","violet")
matplot(times[ota],ero11[ota,],col=colos[1],add=TRUE,type="l",lty=1)
matplot(times[ota],ero22[ota,],col=colos[2],add=TRUE,type="l",lty=1)
matplot(times[ota],ero33[ota,],col=colos[3],add=TRUE,type="l",lty=1)
matplot(times[ota],ero44[ota,],col=colos[4],add=TRUE,type="l",lty=1)
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)
segments(1,0,30000,0,col="green")

Comparison between GARCH and Kernel Regression: Prediction Horizon of One Day

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)
  
file<-"http://jussiklemela.com/statfina/sigmat"
sigmat<-as.matrix(read.csv(file=file))

hoot<-c(0.1,0.2,0.3,0.4,0.5,0.6,0.7)
g1<-20
g2<-10
eta<-1

cspe1<-matrix(0,n,length(hoot))
spe1<-matrix(0,n,length(hoot))
cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

for (hh in 1:length(hoot)){
  h<-hoot[hh]
  #
  x00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     x00[i,1]<-sigmat[i]  
     x00[i,2]<-ma(now,h=g2)
  }
  x000<-copula.trans(x00)
  #
  z00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     z00[i,1]<-sqrt(ma(now^2,h=g1))
     z00[i,2]<-ma(now,h=g2)
  }
  z000<-copula.trans(z00)
  #
  t0<-250
  for (t in (t0+1):(n-eta)){
        hava<-r[t+eta]
        y<-matrix(r[2:t]^2,t-1,1)
        #
        x0<-x000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
        spe1[t,hh]<-abs(sigma^2-hava^2)^2
        cspe1[t,hh]<-sum(spe1[1:t,hh],na.rm=TRUE)
        #
        x0<-z000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
        spe2[t,hh]<-abs(sigma^2-hava^2)^2
        cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)		       
  }
}

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

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

cspe11<-matrix(0,n,1)
minhoot11<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe1[i-1,])
    cspe11[i]<-cspe11[i-1]+spe1[i,ine]
    minhoot11[i]<-hoot[ine]
}
ero11<-cspe11-gcspe[,1]

cspe22<-matrix(0,n,1)
minhoot22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    minhoot22[i]<-hoot[ine]
}
ero22<-cspe22-gcspe[,1]

ota<-c(1000:(n-10))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(2000,-2.0e-05,expression(paste(X[1],"=ewma")),cex=cex,col="orange")
text(2005,-0.0e-05,expression(paste(X[1],"=garch")),cex=cex,col="black")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

plot(times[ota],minhoot11[ota],type="l",lty=1,xlab="",ylab="",col="black")
d<-2
nee<-seq(1,n)  
hee<-(4/(d+2))^(1/(d+4))*nee^(-1/(d+4))
matplot(times[ota],hee[ota],col="blue",add=TRUE,type="l")
mtext("h",side=2,cex=1.5,line=2.4)

ota<-c(1000:(n-6800))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

ota<-c(10000:(n-2000))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(2000,-1.5e-04,"ewma",cex=cex,col="orange")
text(2000,-1.75e-04,"garch",cex=cex,col="black")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

Comparison between GARCH and Kernel Regression: Prediction Horizon of Ten Days

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)
  
file<-"http://jussiklemela.com/statfina/sigmat"
sigmat<-as.matrix(read.csv(file=file))

hoot<-c(0.3,0.4,0.5,0.6,0.7,0.8,0.9)
g1<-20
g2<-20
eta<-10
	   
cspe1<-matrix(0,n,length(hoot))
spe1<-matrix(0,n,length(hoot))
cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

for (hh in 1:length(hoot)){
  h<-hoot[hh]
  #
  x00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     x00[i,1]<-sigmat[i]  
     x00[i,2]<-ma(now,h=g2)
  }
  x000<-copula.trans(x00)
  #
  z00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     z00[i,1]<-sqrt(ma(now^2,h=g1))
     z00[i,2]<-ma(now,h=g2)
  }
  z000<-copula.trans(z00)
  #
  t0<-250
  for (t in (t0+1):(n-eta)){
        hava<-r[t+eta]
        y<-matrix(r[2:t]^2,t-1,1)
        #
        x0<-x000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
        spe1[t,hh]<-abs(sigma^2-hava^2)^2
        cspe1[t,hh]<-sum(spe1[1:t,hh],na.rm=TRUE)
        #
        x0<-z000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
        spe2[t,hh]<-abs(sigma^2-hava^2)^2
        cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)		       
  }
}

# GARCH #########
file<-"http://jussiklemela.com/statfina/garch-netret.csv"
para<-as.matrix(read.csv(file=file))
sigmat<-para[,1]
alpha0t<-para[,2]
alpha1t<-para[,3]		   
betat<-para[,4]
#
gcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,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)
    vest<-barsigma2+(a1+b)^(eta-1)*(sigma^2-barsigma2)
    hava<-r[i+eta]
    diffet[i]<-abs(vest-hava^2)^2
    gcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}

cspe11<-matrix(0,n,1)
minhoot11<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe1[i-1,])
    cspe11[i]<-cspe11[i-1]+spe1[i,ine]
    minhoot11[i]<-hoot[ine]
}
ero11<-cspe11-gcspe

cspe22<-matrix(0,n,1)
minhoot22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    minhoot22[i]<-hoot[ine]
}
ero22<-cspe22-gcspe

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

ota<-c(1000:(n-10))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(2008,-3.0e-05,"ewma",cex=cex,col="orange")
text(2008,-5.5e-05,"garch",cex=cex,col="black")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

plot(times[ota],minhoot11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(0.2,0.95),
cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
d<-2
nee<-seq(1,n)  
hee<-(4/(d+2))^(1/(d+4))*nee^(-1/(d+4))
matplot(times[ota],hee[ota],col="blue",add=TRUE,type="l")
mtext("h",side=2,cex=1.5,line=2.4)

ota<-c(1000:(n-6800))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(2000,-1.5e-04,"ewma",cex=cex,col="orange")
text(2000,-1.75e-04,"garch",cex=cex,col="black")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

ota<-c(10000:(n-2000))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(2000,-1.5e-04,"ewma",cex=cex,col="orange")
text(2000,-1.75e-04,"garch",cex=cex,col="black")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

Comparison between GARCH and Kernel Regression: Return Horizon of Ten Days and Prediction Horizon of One Day

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)

dendat<-r+1
step<-10
leno<-floor(n/step)
data<-matrix(0,leno,1)
for (i in 1:leno){
    cur<-i*step
    sequ<-dendat[(cur-step+1):cur]
    data[i]<-prod(sequ)
}   

r<-data-1
n<-length(r)

# GARCH estimation ########
library("tseries")   
t0<-250
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 (t in (t0+1):length(sigmat)){
    now<-r[1: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)
}

######################

hoot<-c(0.1,0.2,0.3,0.4,0.5,0.6,0.7)
g1<-20
g2<-10
eta<-1

z00<-matrix(0,n,2)
for (i in 1:n){ 
    now<-r[1:i]
    z00[i,1]<-sqrt(ma(now^2,h=g1))
    z00[i,2]<-ma(now,h=g2)
}

x00<-matrix(0,n,2)
for (i in 1:n){ 
    now<-r[1:i]
    x00[i,1]<-sigmat[i]
    x00[i,2]<-ma(now,h=g2)
}
 
cspe1<-matrix(0,n,length(hoot))
spe1<-matrix(0,n,length(hoot))
cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

for (hh in 1:length(hoot)){
  h<-hoot[hh]
  t0<-250
  for (t in (t0+1):(n-eta)){
      z11<-z00[1:t,]
      z000<-copula.trans(z11)
      #
      hava<-r[t+eta]
      y<-matrix(r[2:t]^2,t-1,1)
      #
      x0<-z000[1:t,]
      x<-x0[1:(t-1),]
      arg<-x0[t,]
      est<-max(0,kernesti.regr(arg,x,y,h=h))
      sigma<-sqrt(est)
      spe2[t,hh]<-abs(sigma^2-hava^2)^2
      cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)
      #
      x11<-x00[1:t,]
      x000<-copula.trans(x11)
      #
      x0<-x000[1:t,]
      x<-x0[1:(t-1),]
      arg<-x0[t,]
      est<-max(0,kernesti.regr(arg,x,y,h=h))
      sigma<-sqrt(est)
      spe1[t,hh]<-abs(sigma^2-hava^2)^2
      cspe1[t,hh]<-sum(spe1[1:t,hh],na.rm=TRUE)
  }
}

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

#####################

cspe22<-matrix(0,n,1)
minhoot22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    minhoot22[i]<-hoot[ine]
}
ero22<-cspe22-gcspe
			 
cspe11<-matrix(0,n,1)
minhoot11<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe1[i-1,])
    cspe11[i]<-cspe11[i-1]+spe1[i,ine]
    minhoot11[i]<-hoot[ine]
}
ero11<-cspe11-gcspe

######################
		   
times<-matrix(0,n,1)
delta<-10/251.4
alku<-1950+3*delta
for (i in 1:length(times)) times[i]<-alku+(i-1)*delta

ota<-c(251:(n-10))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(1990,-4.7e-05,expression(paste(X[1],"=ewma")),col="orange")
text(1990,-10.0e-05,expression(paste(X[1],"=garch")),col="black")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

plot(times[ota],minhoot11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(0.2,0.7),
cex=cex,cex.axis=cex.axis,cex.lab=cex.lab)
d<-2
nee<-seq(1,n)  
hee<-(4/(d+2))^(1/(d+4))*nee^(-1/(d+4))
matplot(times[ota],hee[ota],col="blue",add=TRUE,type="l")
mtext("h",side=2,cex=1.5,line=2.4)

Testing statistical significance

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)

############################################
eta<-1

file<-"http://jussiklemela.com/statfina/sigmat"
sigmat<-as.matrix(read.csv(file=file))

hoot<-c(0.1,0.2,0.3,0.4,0.5,0.6,0.7)
g1<-20
g2<-10

cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

for (hh in 1:length(hoot)){
  h<-hoot[hh]
  #
  z00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     z00[i,1]<-sqrt(ma(now^2,h=g1))
     z00[i,2]<-ma(now,h=g2)
  }
  z000<-copula.trans(z00)
  #
  t0<-250
  for (t in (t0+1):(n-eta)){
        hava<-r[t+eta]
        y<-matrix(r[2:t]^2,t-1,1)
        #
        x0<-z000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
        spe2[t,hh]<-abs(sigma^2-hava^2)^2
        cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)		       
  }
}

# garch ########
speg<-matrix(0,n,1)
gcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,1)
for (i in (t0+1):(n-eta)){
     sigma<-sigmat[i]
     hava<-r[i+eta]
     vest<-sigma^2
     diffet[i]<-abs(vest-hava^2)^2
     speg[i]<-abs(vest-hava^2)^2
     gcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}
		   
##
	
cspe22<-matrix(0,n,1)
spe22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    spe22[i]<-spe2[i,ine]
}
ero22<-cspe22-gcspe
eros22<-spe22-speg

##

dee<-ero22
eee<-eros22
      
stepo<-250
ota<-seq(1,n,stepo)
mm<-length(ota)
pute<-0
pets<-matrix(pute,mm,mm)
mins<-10

for (i in 1:(mm-1)){
for (j in (i+1):mm){
    t00<-ota[i]
    len<-ota[j]
    otas<-seq(t00,len)
    sigts<-var(eee[t00:len],na.rm=TRUE)
    norma<-(len-t00+1)^(-1/2)
    zum<-dee[len]-dee[t00]
    ts11<-norma*zum/sqrt(sigts)
    pets[i,j]<-pnorm(ts11)
    mins<-min(mins,pets[i,j])
}}

###################################

eta<-10

hoot<-c(0.3,0.4,0.5,0.6,0.7,0.8,0.9)
g1<-20
g2<-20
eta<-10
	   
cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

for (hh in 1:length(hoot)){
  h<-hoot[hh]
  #
  z00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     z00[i,1]<-sqrt(ma(now^2,h=g1))
     z00[i,2]<-ma(now,h=g2)
  }
  z000<-copula.trans(z00)
  #
  t0<-250
  for (t in (t0+1):(n-eta)){
        hava<-r[t+eta]
        y<-matrix(r[2:t]^2,t-1,1)
        #
        x0<-z000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
        spe2[t,hh]<-abs(sigma^2-hava^2)^2
        cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)		       
  }
}

# GARCH #########
file<-"http://jussiklemela.com/statfina/garch-netret.csv"
para<-as.matrix(read.csv(file=file))
sigmat<-para[,1]
alpha0t<-para[,2]
alpha1t<-para[,3]		   
betat<-para[,4]
#
speg<-matrix(0,n,1)
gcspe<-matrix(0,n,1)
diffet<-matrix(0,n-eta,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)
    vest<-barsigma2+(a1+b)^(eta-1)*(sigma^2-barsigma2)
    hava<-r[i+eta]
    diffet[i]<-abs(vest-hava^2)^2
    speg[i]<-abs(vest-hava^2)^2	
    gcspe[i]<-sum(diffet[1:i],na.rm=TRUE)
}

##

cspe22<-matrix(0,n,1)
spe22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    spe22[i]<-spe2[i,ine]
}
ero22<-cspe22-gcspe[,1]
eros22<-spe22-speg

##

dee<-ero22
eee<-eros22
      
stepo<-250
ota<-seq(1,n,stepo)
mm<-length(ota)
pute<-0
pets10<-matrix(pute,mm,mm)
mins<-10

for (i in 1:(mm-1)){
for (j in (i+1):mm){
    t00<-ota[i]
    len<-ota[j]
    otas<-seq(t00,len)
    sigts<-var(eee[t00:len],na.rm=TRUE)
    norma<-(len-t00+1)^(-1/2)
    zum<-dee[len]-dee[t00]
    ts11<-norma*zum/sqrt(sigts)
    pets10[i,j]<-pnorm(ts11)
    mins<-min(mins,pets[i,j])
}}

#######################################
	     
alku<-1
times<-matrix(0,mm,1)
delta<-1/252
alko<-1953+5*delta+alku*delta  #stepo*delta
for (i in 1:length(times)) times[i]<-alko+(i-1)*delta*stepo

good<-matrix(2,mm,mm)      
hyva<-which(pets<=0.05,arr.ind=TRUE)
good[hyva]<-1
kesk<-which(pets<=0.01,arr.ind=TRUE)
good[kesk]<-0		    
alat<-which(pets<=0.0,arr.ind=TRUE)
good[alat]<--1	    
image(times,times,good,col=c("white","red","orange","blue"),
xlab="",ylab="")
mtext("end of the period",side=2,line=2.4)
mtext("beginning of the period",side=1,line=2.4)

good<-matrix(2,mm,mm)      
hyva<-which(pets10<=0.05,arr.ind=TRUE)
good[hyva]<-1
kesk<-which(pets10<=0.01,arr.ind=TRUE)
good[kesk]<-0		    
alat<-which(pets10<=0.0,arr.ind=TRUE)
good[alat]<--1	    
image(times,times,good,col=c("white","red","orange","blue"),
xlab="",ylab="")
mtext("end of the period",side=2,line=2.4)
mtext("beginning of the period",side=1,line=2.4)
  

Contour plot and slices

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
times0<-matrix(0,length(y),1)
delta<-1/251.63
alku<-1950+2/12
for (i in 1:length(times0)) times0[i]<-alku+(i-1)*delta

m<-length(y)	
r<-y[2:m]/y[1:(m-1)]-1
times<-times0[2:m]
n<-length(r)

h<-0.5
eta<-1

h2<-40
x00<-matrix(0,n,2)
for (i in 1:n){ 
   now<-r[1:i]
   x00[i,1]<-sqrt(ma(now^2,h=h2))
   x00[i,2]<-ma(now,h=h2)
}
x0<-copula.trans(x00)

y<-matrix(r[2:n]^2,n-1,1)
x<-x0[1:(n-1),]
min1<-min(x[,1])
max1<-max(x[,1])
min2<-min(x[,2])
max2<-max(x[,2])

num<-100
step1<-(max1-min1)/num
x1<-seq(min1,max1,step1)
step2<-(max2-min2)/num
x2<-seq(min2,max2,step2)
z<-matrix(0,length(x1),length(x2))
for (i in 1:length(x1)){
   for (j in 1:length(x2)){
      arg<-c(x1[i],x2[j])
      z[i,j]<-max(0,kernesti.regr(arg,x,y,h=h))
   }
}

num<-50
step1<-(max1-min1)/num
x11<-seq(min1,max1,step1)
step2<-(max2-min2)/num
x22<-seq(min2,max2,step2)
z2<-matrix(0,length(x11),length(x22))
for (i in 1:length(x11)){
   for (j in 1:length(x22)){
      arg<-c(x11[i],x22[j])
      z2[i,j]<-max(0,kernesti.regr(arg,x,y,h=h))
   }
}

plot(x,col="yellow",xlab="",ylab="")
contour(x1,x2,z,nlevels=20,add=TRUE)
mtext("MA of squared returns",side=1,line=2.4)
mtext("MA of returns",side=2,line=2.4)

persp(x11,x22,z2,theta=-30,phi=20,ticktype="detailed",
xlab="MA of squared returns",ylab="MA of returns",zlab="") 

plot(x1,z[,1],type="l",xlab="",ylab="")
for (i in 1:length(x2)){
    matplot(x1,z[,i],add=TRUE,type="l")
}
mtext("MA of squared returns",side=1,line=2.4)

plot(x22,z2[1,],type="l",xlab="",ylab="",ylim=c(min(z2),max(z2)))
for (i in 1:length(x11)){
    matplot(x11,z2[i,],add=TRUE,type="l")
}
mtext("MA of returns",side=1,line=2.4)

xeka<-matrix(0,length(x1),1)
deda<-x00[,1]
for (i in 1:length(xeka)){	     
    pee<-pnorm(x1)[i]
    quu<-sort(deda)[round(pee*length(deda))]
    xeka[i]<-quu
}
plot(xeka,z[,1],type="l",xlab="",ylab="")
for (i in 1:length(x2)){
    matplot(xeka,z[,i],add=TRUE,type="l")
}
mtext("MA of squared returns",side=1,line=2.4)

xtok<-matrix(0,length(x22),1)
deda<-x00[,2]
for (i in 1:length(xtok)){	     
    pee<-pnorm(x22)[i]
    quu<-sort(deda)[round(pee*length(deda))]
    xtok[i]<-quu
}
plot(xtok,z2[1,],type="l",xlab="",ylab="",ylim=c(min(z2),max(z2)))
for (i in 1:length(xtok)){
    matplot(xtok,z2[i,],add=TRUE,type="l")
}
mtext("MA of returns",side=1,line=2.4)

oni<-10
xtok[oni]
plot(xtok,z2[oni,],type="l",xlab="",ylab="")
mtext("MA of returns",side=1,cex=cex,line=2.4)
		

Quantile Estimation with One Day Horizon

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)
  
rho<-function(x,p)
{
if (x<0) ans<-x*(p-1) else ans<-x*p
return(ans)
}    
    
hoot<-c(0.1,0.2,0.3,0.4,0.5,0.6,0.7)
g1<-20
g2<-10
eta<-1
alpha<-0.01
	      
cspe1<-matrix(0,n,length(hoot))
spe1<-matrix(0,n,length(hoot))
cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

sigmat.oma<-matrix(0,n,1)	

for (hh in 1:length(hoot)){
  h<-hoot[hh]
  #
  z00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     z00[i,1]<-sqrt(ma(now^2,h=g1))
     z00[i,2]<-ma(now,h=g2)
  }
  z000<-copula.trans(z00)
  #
  t0<-250
  for (t in 1:t0){
        now<-r[1:t]
        hava<-r[t+eta]
        if (t>1) sigma<-sd(now) else sigma<-1
        sigmat.oma[t]<-sigma
  } 
  for (t in (t0+1):(n-eta)){
        now<-r[1:t]		      
        hava<-r[t+eta]
        y<-matrix(r[2:t]^2,t-1,1)
        #
        x0<-z000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
	sigmat.oma[t]<-sigma
        #
	df<-5       
        ker<-qt(alpha,df=df)*sqrt((df-2)/df)
        mui<-mean(now)
        qhat<-mui+sigma*ker
        spe1[t,hh]<-rho(hava-qhat,alpha)  #abs(sigma^2-hava^2)^2
        cspe1[t,hh]<-sum(spe1[1:t,hh],na.rm=TRUE)
        #
        residu<-now[t0:t]/sigmat.oma[t0:t]
        ncur<-length(residu)
        m<-ceiling(ncur*alpha)
        eqvan<-sort(residu)[m]
        qhat<-mui+sigma*eqvan
        spe2[t,hh]<-rho(hava-qhat,alpha)  
	cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)
   }
}

# garch #########
file<-"http://jussiklemela.com/statfina/sigmat"
sigmat<-as.matrix(read.csv(file=file))
#
cspe01<-matrix(0,n,1)
spe01<-matrix(0,n,1)
cspe02<-matrix(0,n,1)
spe02<-matrix(0,n,1)
eta<-1
for (t in (t0+1):(n-eta)){
    now<-r[1:t]
    hava<-r[t+eta]
    sigma<-sigmat[t]
    #			      
    df<-5       
    ker<-qt(alpha,df=df)*sqrt((df-2)/df)
    mui<-mean(now)
    qhat<-mui+sigma*ker
    spe01[t]<-rho(hava-qhat,alpha)  #abs(sigma^2-hava^2)^2
    cspe01[t]<-sum(spe01[1:t],na.rm=TRUE)
    #
    residu<-now[t0:t]/sigmat[t0:t]
    ncur<-length(residu)
    m<-ceiling(ncur*alpha)
    eqvan<-sort(residu)[m]
    qhat<-mui+sigma*eqvan
    spe02[t]<-rho(hava-qhat,alpha)  
    cspe02[t]<-sum(spe02[1:t],na.rm=TRUE)
}

#####################

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

###################

cspe11<-matrix(0,n,1)
minhoot11<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe1[i-1,])
    cspe11[i]<-cspe11[i-1]+spe1[i,ine]
    minhoot11[i]<-hoot[ine]
}
ero11<-cspe11-cspe01

cspe22<-matrix(0,n,1)
minhoot22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    minhoot22[i]<-hoot[ine]
}
ero22<-cspe22-cspe02

ota<-c(1000:(n-10))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(1970,-0.01,"student",col="black")
text(1970,-0.10,"emp",col="orange")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

plot(times[ota],minhoot11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(0.2,0.7))
matplot(times[ota],minhoot22[ota],type="l",lty=1,add=TRUE,col="orange")     
d<-2
nee<-seq(1,n)  
hee<-(4/(d+2))^(1/(d+4))*nee^(-1/(d+4))
matplot(times[ota],hee[ota],col="blue",add=TRUE,type="l")
mtext("h",side=2,cex=1.5,line=2.4)

Quantile Estimation with Ten Day Returns

file<-"http://jussiklemela.com/statfina/sp500"
y<-read.csv(file=file)[,1]
m<-length(y)
r<-(y[2:m]-y[1:(m-1)])/y[1:(m-1)]  #r<-log(y[2:m])-log(y[1:(m-1)])
n<-length(r)

dendat<-r+1
step<-10
leno<-floor(n/step)
data<-matrix(0,leno,1)
for (i in 1:leno){
    cur<-i*step
    sequ<-dendat[(cur-step+1):cur]
    data[i]<-prod(sequ)
}   
r<-data-1
n<-length(r)

rho<-function(x,p)
{
if (x<0) ans<-x*(p-1) else ans<-x*p
return(ans)
}    
    
hoot<-seq(0.1,1.8,0.1)
g1<-20
g2<-10
eta<-1
alpha<-0.01
	      
cspe1<-matrix(0,n,length(hoot))
spe1<-matrix(0,n,length(hoot))
cspe2<-matrix(0,n,length(hoot))
spe2<-matrix(0,n,length(hoot))

sigmat.oma<-matrix(0,n,1)
t0<-250
		  
for (hh in 1:length(hoot)){
  h<-hoot[hh]
  #
  z00<-matrix(0,n,2)
  for (i in 1:n){ 
     now<-r[1:i]
     z00[i,1]<-sqrt(ma(now^2,h=g1))
     z00[i,2]<-ma(now,h=g2)
  }
  z000<-copula.trans(z00)
  #
  for (t in (t0+1):(n-eta)){
        now<-r[1:t]		      
        hava<-r[t+eta]
        y<-matrix(r[2:t]^2,t-1,1)
        #
        x0<-z000[1:t,]
        x<-x0[1:(t-1),]
        arg<-x0[t,]
        est<-max(0,kernesti.regr(arg,x,y,h=h))
        sigma<-sqrt(est)
	sigmat.oma[t]<-sigma
        #
	df<-7 
        ker<-qt(alpha,df=df)*sqrt((df-2)/df)
        mui<-mean(now)
        qhat<-mui+sigma*ker
        spe1[t,hh]<-rho(hava-qhat,alpha)  #abs(sigma^2-hava^2)^2
        cspe1[t,hh]<-sum(spe1[1:t,hh],na.rm=TRUE)
        #
        residu<-now[t0:t]/sigmat.oma[t0:t]
        ncur<-length(residu)
        m<-ceiling(ncur*alpha)
        eqvan<-sort(residu)[m]
        qhat<-mui+sigma*eqvan
        spe2[t,hh]<-rho(hava-qhat,alpha)  
	cspe2[t,hh]<-sum(spe2[1:t,hh],na.rm=TRUE)
   }
}

# GARCH estimation ########
library("tseries")   
t0<-250
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 (t in (t0+1):length(sigmat)){
    now<-r[1: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)
}
#####
cspe01<-matrix(0,n,1)
spe01<-matrix(0,n,1)
cspe02<-matrix(0,n,1)
spe02<-matrix(0,n,1)
for (t in (t0+1):(n-eta)){
    now<-r[1:t]
    hava<-r[t+eta]
    sigma<-sigmat[t]
    #			            
    ker<-qt(alpha,df=df)*sqrt((df-2)/df)
    mui<-mean(now)
    qhat<-mui+sigma*ker
    spe01[t]<-rho(hava-qhat,alpha)  #abs(sigma^2-hava^2)^2
    cspe01[t]<-sum(spe01[1:t],na.rm=TRUE)
    #
    residu<-now[t0:t]/sigmat[t0:t]
    ncur<-length(residu)
    m<-ceiling(ncur*alpha)
    eqvan<-sort(residu)[m]
    qhat<-mui+sigma*eqvan
    spe02[t]<-rho(hava-qhat,alpha)  
    cspe02[t]<-sum(spe02[1:t],na.rm=TRUE)
}

###################

cspe11<-matrix(0,n,1)
minhoot11<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe1[i-1,])
    cspe11[i]<-cspe11[i-1]+spe1[i,ine]
    minhoot11[i]<-hoot[ine]
}
ero11<-cspe11-cspe01

cspe22<-matrix(0,n,1)
minhoot22<-matrix(0,n,1)
for (i in 2:n){
    ine<-which.min(cspe2[i-1,])
    cspe22[i]<-cspe22[i-1]+spe2[i,ine]
    minhoot22[i]<-hoot[ine]
}
ero22<-cspe22-cspe02

#####################

times<-matrix(0,n,1)
delta<-10/251.4
alku<-1950+3*delta
for (i in 1:length(times)) times[i]<-alku+(i-1)*delta
	  
ota<-c(250:(n-10))
matplot(times[ota],ero11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(min(ero11[ota],ero22[ota]),max(ero11[ota],ero22[ota])))
matplot(times[ota],ero22[ota],type="l",lty=1,add=TRUE,col="orange")
segments(0,0,20000,0,col="green")
text(1970,0.04,"student",col="black")
text(1970,-0.05,"emp",col="orange")
mtext("differenced CSPE",side=2,cex=1.5,line=2.4)

plot(times[ota],minhoot11[ota],type="l",lty=1,xlab="",ylab="",col="black",
ylim=c(0.2,1.9))
matplot(times[ota],minhoot22[ota],type="l",lty=1,add=TRUE,col="orange")     
d<-2
nee<-seq(1,n)  
hee<-(4/(d+2))^(1/(d+4))*nee^(-1/(d+4))
matplot(times[ota],hee[ota],col="blue",add=TRUE,type="l")
mtext("h",side=2,cex=1.5,line=2.4)
	   

Jussi Klemelä

July 2017