R code to verify Volatility Forecast Theory by Robert Engle

The Paper from Christian Brownlees, Rob Engle and Bryan Kelly can be downloaded from here:

http://faculty.chicagobooth.edu/bryan.kelly/research/pdf/volfor.pdf‎

Based on the real market data of SPX, I wrote the following R code to verify that GARCH model ranking doesn't change as the forecast horizon changes:

rm(list=ls())

library(rugarch)

library(timeSeries)

library(xts)

fname="spx.csv";

res=readSeries(fname,header=T,sep=",",format="%m/%d/%Y");

colnames(res)=c("O", "H", "L", "C", "V", "A")

sp500=as.xts(res[,6]);

rtn=na.omit(returns(sp500));

rtnsquare=rtn^2;

ord=ar(as.numeric(rtn),method='mle');

armaorder_=c(ord$order,0)

spec1=ugarchspec(mean.model = list(armaOrder = armaorder_),

                 variance.model = list(model='fGARCH', submodel='GARCH'));

spec2=ugarchspec(mean.model = list(armaOrder = armaorder_),

                 variance.model = list(model='fGARCH', submodel='NGARCH'));

spec3=ugarchspec(mean.model = list(armaOrder = armaorder_),

                 variance.model = list(model='fGARCH', submodel='TGARCH'));

spec4=ugarchspec(mean.model = list(armaOrder = armaorder_),

                 variance.model = list(model='eGARCH', submodel=NULL));

spec5=ugarchspec(mean.model = list(armaOrder = armaorder_),

                 variance.model = list(model='fGARCH', submodel='APARCH'));

#Quasi Loss Function

QL=function(squaredR,VarForecast){

  tmp=squaredR/VarForecast;

  return(tmp-log(tmp)-1);

}

specs=c(spec1,spec2,spec3,spec4,spec5)

sname=c("GARCH","NGARCH","TARCH","eGARCH","APARCH")

Cal=function(specs){

  horizons=c(1,10,20)

  for(horizon in horizons){

    tmp=0

    count=0;

    index=1;

    for(spec_ in specs){

      cat(sname[index],": ")

      index=index+1;

      out_sample_len=c(50:50);

      for(i in out_sample_len){

        fit1=ugarchfit(rtn, spec = spec_, out.sample=i,solver='hybrid')

        fcast1=ugarchforecast(fit1,n.ahead=horizon)

        squaredR=rtnsquare[dim(rtnsquare)[1]-i+horizon];

        VarForecast=sigma(fcast1)[horizon]^2;

        loss=QL(squaredR,VarForecast);

        tmp=tmp+loss;

        count=count+1;

      }

      cat("Horizon =",horizon,"; LQ =",(tmp/count),"\n");

    }

  }

}

Cal(specs)

The result is as follows:

GARCH : Horizon = 1 ; LQ = 1.234608
NGARCH : Horizon = 1 ; LQ = 1.091949
TARCH : Horizon = 1 ; LQ = 1.067449
eGARCH : Horizon = 1 ; LQ = 1.060457
APARCH : Horizon = 1 ; LQ = 1.02967

GARCH : Horizon = 10 ; LQ = 2.259468
TARCH : Horizon = 10 ; LQ = 2.123502
NGARCH : Horizon = 10 ; LQ = 2.141799
eGARCH : Horizon = 10 ; LQ = 2.107577
APARCH : Horizon = 10 ; LQ = 2.089418

GARCH : Horizon = 20 ; LQ = 1.25515
NGARCH : Horizon = 20 ; LQ = 1.179098
TARCH : Horizon = 20 ; LQ = 1.170422
eGARCH : Horizon = 20 ; LQ = 1.150264
APARCH : Horizon = 20 ; LQ = 1.144654

The GarchFit takes very long time to finish as my notebook is not so powerful. I might need a cluster or grid to do this kind of verification.