vech <- function (x)
{
t(x[!upper.tri(x)])
}
fm2 = function(u){
return( mean(u^2))
}
fm3 = function(u){
return( mean(u^3))
}
fm4 = function(u){
return( mean(u^4))
}
fm6 = function(u){
return( mean(u^6))
}
# U = matrix( rnorm(400) , ncol = 4 )
coskewCC = function( U ){
# Implementation for N >= 3
# The coskewness matrix is based on stacking subcomoment matrices columnwise
# Advantage: matrix notation:
# M3 = E[ UU' %x% U ]: [U[1]%*%UU' | U[2]%*%UU' | ... | U[N]%*%UU' ]
# element [j,(i-1)N+k] is E[U[i]*U[j]*U[k]]
# U needs to be a TxN matrix holding the centred returns of asset i in column i
# N rows, N^3 columns
# See Martellini and Ziemann (Review of Financial Studies)
N = ncol(U); T = nrow(U);
vm2 = apply(U , 2 , 'fm2' )
vm3 = apply(U , 2 , 'fm3' )
vm4 = apply(U , 2 , 'fm4' )
vr2 = c();
for( j in 2:N ){
for( i in 1:(j-1) ){
vr2 = c( vr2 , sum( U[,i]^2*U[,j] )/sqrt( vm4[i]*vm2[j] ) )
}
}
r2 = sum(vr2)*2/( T*N*(N-1) )
vr5 = c();
for( j in 2:N ){
for( i in 1:(j-1) ){
vr5 = c( vr5 , sum( U[,i]^2*U[,j]^2 )/sqrt( vm4[i]*vm4[j] ) )
}
}
r5 = sum(vr5)*2/( T*N*(N-1) )
vr4 = c();
for( k in 3:N ){
for( j in 2:(k-1) ){
for( i in 1:(j-1) ){
#vr4 = c( vr4 , sum( U[,i]*U[,j]*U[,k] )/sqrt( vm2[k]*r5*sqrt(vm4[i]*vm4[j]) ) )
denom = mean(
c(vm2[k]*r5*sqrt( vm4[i]*vm4[j] ) , vm2[i]*r5*sqrt( vm4[j]*vm4[k] ) , vm2[j]*r5*sqrt( vm4[i]*vm4[k] ) ) )
vr4 = c( vr4 , sum( U[,i]*U[,j]*U[,k] )/sqrt( denom ) )
}
}
}
r4 = sum(vr4)*6/( T*N*(N-1)*(N-2) )
M3CC = matrix( rep(0,N^3) , nrow = N )
for( i in 1:N ){
for( j in 1:N ){
for( k in 1:N ){
if( (i==j) & (j==k) & (i==k) ){
sijk = vm3[i]
}else{
if( (i==j) | (i==k) | (j==k) ){
if( (i==j) ){ sijk = r2*sqrt( vm4[i]*vm2[k] )}
if( (i==k) ){ sijk = r2*sqrt( vm4[i]*vm2[j] )}
if( (j==k) ){ sijk = r2*sqrt( vm4[j]*vm2[i] )}
}else{
# sijk = r4*sqrt( vm2[k]*r5*sqrt( vm4[i]*vm4[j] ) )
sijk_all = r4*sqrt(
mean(
c(vm2[k]*r5*sqrt( vm4[i]*vm4[j] ) , vm2[i]*r5*sqrt( vm4[j]*vm4[k] ) , vm2[j]*r5*sqrt( vm4[i]*vm4[k] ) ) ))
sijk = sijk_all
}
}
# i denotes the subcomoment matrix
M3CC[j,(i-1)*N+k] = sijk;
}
}
}
return(M3CC)
}
# cokurtCC(U)
cokurtCC = function( U ){
# Implementation for N >=4
# The coskewness matrix is based on stacking subcomoment matrices columnwise
# Advantage: matrix notation
# M4 = E[ UU' %x% U %x% U ]:
#[ U[1]^2%*%UU'| U[1]U[2]%*%UU' | .... | U[1]U[N]%*%UU'|
# U[2]U[1]%*%UU'| U[2]^2%*%UU' | .... | U[2]U[N]%*%UU'|
# ...
# U[N]U[1]%*%UU'| U[N]U[2]%*%UU' | .... | U[N]^2%*%UU ]
# element [k,( (i-1)*(N^2)+(j-1)*N+l )] is E[ U[i]*U[j]*U[k]*U[l] ]
# U needs to be a TxN matrix holding the centred returns of asset i in column i
# N rows, N^4 columns
# See Martellini and Ziemann (Review of Financial Studies)
N = ncol(U); T = nrow(U);
vm2 = apply(U , 2 , 'fm2' )
vm3 = apply(U , 2 , 'fm3' )
vm4 = apply(U , 2 , 'fm4' )
vm6 = apply(U , 2 , 'fm6' )
vr3 = c();
for( j in 2:N ){
for( i in 1:(j-1) ){
vr3 = c( vr3 , sum( U[,i]^3*U[,j] )/sqrt( vm6[i]*vm2[j] ) )
}
}
r3 = sum(vr3)*2/( T*N*(N-1) )
vr5 = c();
for( j in 2:N ){
for( i in 1:(j-1) ){
vr5 = c( vr5 , sum( U[,i]^2*U[,j]^2 )/sqrt( vm4[i]*vm4[j] ) )
}
}
r5 = sum(vr5)*2/( T*N*(N-1) )
vr6 = c();
for( k in 3:N ){
for( j in 2:(k-1) ){
for( i in 1:(j-1) ){
vr6 = c( vr6 , sum( U[,i]^2*U[,j]*U[,k] )/sqrt( vm4[i]*r5*sqrt(vm4[j]*vm4[k]) ) )
}
}
}
r6 = sum(vr6)*6/( T*N*(N-1)*(N-2) )
vr7 = c();
for( l in 4:N ){
for( k in 3:(l-1) ){
for( j in 2:(k-1) ){
for( i in 1:(j-1) ){
vr7 = c( vr7 , sum( U[,i]*U[,j]*U[,k]*U[,l] )/sqrt( r5*sqrt(vm4[i]*vm4[j])*r5*sqrt(vm4[k]*vm4[l]) ) )
}
}
}
}
r7 = sum(vr7)*24/( T*N*(N-1)*(N-2)*(N-3) )
M4CC = matrix( rep(0,N^4) , nrow = N )
for( i in 1:N ){
for( j in 1:N ){
for( k in 1:N ){
for( l in 1:N ){
# in the general case: i, j , k and l are all different
kijkl = r7*sqrt( r5*sqrt(vm4[i]*vm4[j])*r5*sqrt(vm4[k]*vm4[l]) )
# if at least one of them are equal:
if( (i==j) | (i==k) | (i==l) | (j==k) | (j==l) | (k==l) ){
if( (i==j) & (i==k) & (i==l) ){
# These are the kurtosis estimates of the individual assets: E[u^4]
kijkl = vm4[i];
}else{
if( ((i==j) & (i==k)) | ((i==j) & (i==l)) | ((i==k) & (i==l)) | ((j==k) & (j==l)) ){
# kiij E[ U[,i]^3*U[,j] ] = r3*sqrt( vm6[i]*vm2[j] )
if( ((i==j) & (i==k)) ){ kijkl = r3*sqrt( vm6[i]*vm2[l] ) }
if( ((i==j) & (i==l)) ){ kijkl = r3*sqrt( vm6[i]*vm2[k] ) }
if( ((i==k) & (i==l)) ){ kijkl = r3*sqrt( vm6[i]*vm2[j] ) }
if( ((j==k) & (j==l)) ){ kijkl = r3*sqrt( vm6[j]*vm2[i] ) }
}else{
if( ((i==j) & (k==l))| ( (i==k) & (j==l) ) | ( (i==l) & (j==k)) ){
# kiijj = E[ U[,i]^2*U[,j]^2 ] = r5*sqrt( vm4[i]*vm4[j] )
if ( ((i==j) & (k==l)) ){ kijkl = r5*sqrt( vm4[i]*vm4[k] ) }
if ( ((i==k) & (j==l)) ){ kijkl = r5*sqrt( vm4[i]*vm4[j] ) }
if ( ((i==l) & (j==k)) ){ kijkl = r5*sqrt( vm4[i]*vm4[j] ) }
}else{
# kiijk = E[ U[,i]^2*U[,j]*U[,k] ] = r6*sqrt( vm4[i]*r5*sqrt( vm4[j]*vm4[k] ) )
if ( (i==j) ){ kijkl = r6*sqrt( vm4[i]*r5*sqrt(vm4[k]*vm4[l]) ) }
if ( (i==k) ){ kijkl = r6*sqrt( vm4[i]*r5*sqrt(vm4[j]*vm4[l]) ) }
if ( (i==l) ){ kijkl = r6*sqrt( vm4[i]*r5*sqrt(vm4[j]*vm4[k]) ) }
if ( (j==k) ){ kijkl = r6*sqrt( vm4[j]*r5*sqrt(vm4[i]*vm4[l]) ) }
if ( (j==l) ){ kijkl = r6*sqrt( vm4[j]*r5*sqrt(vm4[i]*vm4[k]) ) }
if ( (k==l) ){ kijkl = r6*sqrt( vm4[k]*r5*sqrt(vm4[i]*vm4[j]) ) }
}
}
}
}
M4CC[k,( (i-1)*(N^2)+(j-1)*N+l )] = kijkl;
}#loop l
}#loop k
}#loop j
}#loop i
return(M4CC)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.