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)