Feeds:
Posts

## KURTOSIS OF SUBORDINATION BY POWER-LAW WAITING TIMES

When Peter Clark wrote his famous paper on explaining kurtosis in financial returns by Bochner’s subordination in 1971, he did not have an exact model for waiting times.  We know today that these waiting times have power law tails and the machinery built for anomalous diffusions gives us an exact model depending on the power law parameter.  So how good are the predicted kurtosis of returns by the following simple model: calculate the tail index of the waiting times directly for a stock and then use the subordination formula on a Gaussian kernel to obtain a subordinated kernel for returns, i.e. Brownian motion subordinated to an inverse stable subordinator (see meerschaert-inverse-stable-subordinator-2011).  Fix time $t=1.0$ and just calculate the kurtosis of the distribution thus obtained.  So this is a direct model for kurtosis.  As a test run for this quantitative kurtosis prediction model, we define outlier as any return value greater than 8 standard deviations so this is a limited check.  What are the relative errors of this method?

Over 42 stocks the outlier-adjusted kurtosis have mean $8.55 \pm 3.36$ while the predictions (by method described above) have mean $9.85 \pm 0.2$ with relative errors $0.47 \pm 0.5$. Now this is far from a perfect model but given the simplicity of our model — we are predicting kurtosis based only on waiting time power law for volatility jumps — this is actually reasonable.

Stock,OutlierAdjustedKurtosis,PredictedKurtosis,RelErr
LRN,9.83141668515334,10.0605389835161,0.0233051151934965
PRE,7.50405788304994,10.0904904144613,0.344671186139639
IRBT,10.4320581565518,9.67739760058358,0.0723405242420177
G,6.44990465506245,9.97469350160842,0.546486969195645
ICUI,7.27998231715923,9.69960985025275,0.332367226688223
IGT,5.57760494842199,9.72694116789355,0.743927950767735
SIAL,7.97547383524259,9.88458731734433,0.239373048114785
KND,12.3165475183667,9.96703599044731,0.190760562114974
ASH,11.0928719361077,9.87656481375115,0.109647630420886
ABMD,5.46714688189871,9.5450826347754,0.74589833435396
WSM,6.98752035950721,9.65067691759533,0.381130418384347
RTEC,5.29297240180037,9.83686940836542,0.858477366142977
LXU,11.9564165703199,10.3624468451538,0.133315004189707
ELNK,8.31628672113459,9.7460337653134,0.171921326443126
CVO,6.28560820543105,9.58309283985844,0.524608681714875
GOV,2.5941242060119,10.1634035673209,2.91785541485143
RRGB,9.60436370848588,10.0276626767439,0.0440736087372495
OEH,12.2364959644304,9.9111301801339,0.190035267535409
GDP,5.94844587927817,9.91158655574775,0.666248085113366
SIRO,7.74410306855586,10.3348226229963,0.334540944445814
FFIN,8.50812557408157,9.75064710958655,0.14603939783048
PBCT,9.11669260316958,9.56881611736894,0.0495929317658662
SRI,9.49786804031522,9.83184046414371,0.0351628831239694
MM,16.0170924928767,9.81242901171704,0.387377639476022
SWKS,5.66164093002866,9.5123038101636,0.68013194897463
TOL,4.01655476166821,9.61004415858085,1.39260877264611
HBI,19.5742437522558,9.8621176238896,0.496168651585681
MCRI,7.82861359045437,9.87615774095122,0.261546201870721
WABC,6.69867481853173,10.0559447165258,0.501184187759983
ITRI,11.1667563031087,9.73218636625486,0.128467918338514
APH,6.5997736237985,9.83962788899471,0.490903847597655
DF,12.4731221164775,9.97122464445926,0.200583097692371
HTZ,12.4116888716239,10.0290428672707,0.191967912586054
KNL,3.97213972146467,9.94084220850729,1.50264162531567
ACGL,8.57218294419001,9.8569376498777,0.149874858487296
SCSS,11.8279977668077,9.65948549584479,0.183337223570357
MRK,5.46177969794201,10.1016745529448,0.849520689520137
BKD,13.3226501856739,9.89219452482362,0.257490485229369
TITN,6.5678578552141,9.69240534348549,0.475733116816904
FNSR,5.90233914008764,9.56431517305507,0.620427926293822
CR,6.41159193301708,9.82430736983747,0.532272713621448


Full code to generate these results:

library(e1071)
library(stabledist)
library(Quandl)
library(condmixt)
library(igraph)

Quandl.api_key("CoHxdmas1Pz2t7Hk8maE")

cap.outliers<-function(rets){
r<-rets
outlier.threshold<-8*sd(rets)
r[r< -outlier.threshold] outlier.threshold]<- outlier.threshold
return(r)
}

voljump.tail <-
function(vx){
dvx<-diff(vx)
N<-length(dvx)
jumps<-rep(0,N)
threshthresh){
jumps[k]<-1.0
}
}
const<-1
x<-c()
for (k in 1:N){
if (jumps[k]==0.0){
const<-const+1
}else{
x<-append(x,const)
const<-1
}
}
return(x)
}

kurtosis.MLWT <-
function(b,t){
dr=0.5 #changed 11/15/2012 to show peak at zero
r=seq(dr,5000.0,dr) #changed 11/15/2012 to show peak at zero
#b=0.75
pi=3.1415927
g=(cos(pi*b/2))^(1/b)
h=dstable(r, alpha = b, beta = 1.0, gamma = g, delta = 0, pm=1)
D=1.0
mcall <- function(y,t) {
sum(dnorm(y, mean = 0.0, sd =sqrt(D*(t/r)^b) )*h*dr)
}
x=seq(-5.0,5.0,0.1)
m=x
#t=0.1
for (i in 1:length(x)){
m[i]=mcall(x[i],t)}
#plot(x,m,type="l")
integrand2<-function(x){x^2*mcall(x,t)}
integrand4<-function(x){x^4*mcall(x,t)}
m2<-integrate(integrand2,lower=0,upper=Inf,rel.tol=1e-2)
m4<-integrate(integrand4,lower=0,upper=Inf,rel.tol=1e-2)
return(m4$value/m2$value^2)
}

compare.kurtosis <-
function(stock){
d<-Quandl(paste("WIKI/",stock,sep=""))
r<-diff(log(px))
r[is.na(r)]<-0
vx<-log(r^2+1e-6)
jumptimes<-voljump.tail(vx)
actual.kurtosis<-kurtosis(cap.outliers(r))

beta.kur10<-1.0/(1.0*power.law.fit(jumptimes,xmin=median(jumptimes))\$alpha)
#times<-c(100.0)
#times<-c(1.0,5.0,10.0,20.0,40.0,80.0)
#for (t in times){
#    print(paste('t=',t,'k=',kurtosis.MLWT(beta.kur10,t)))
#}
predicted.kurtosis<-kurtosis.MLWT(beta.kur10,20.0)
relative.error<-abs(actual.kurtosis-predicted.kurtosis)/actual.kurtosis
print(paste(stock,'act kur=',actual.kurtosis,'pred=',predicted.kurtosis,' relerr=',relative.error))
#print(paste(' k10=', k10,' k15=',k15,' k20=', k20,' k25=',k25,' k30=',k30,' k35=',k35,' k40=',k40,' k45=',k45,' k50=',k50))

}