DAX-FTSE copula plot

# This script requires the package 'copula' to be installed.

#
# copula DAX, FTSE
#

x = EuStockMarkets[,1]  # DAX
x = diff( log(x) )
y = EuStockMarkets[,4]  # FTSE
y = diff( log(y) )

# helpers for copula stuff: qkdens, pkdens

scdf1 = function( t, x, h ) {
    mean( pnorm( (t-x)/h ) )
}
scdf = Vectorize( scdf1, "t" )

pkdensx1 = function( z, mean=0, sd=1 ) {
  h = bw.bcv( x )
  return( scdf1( z, x, h ) )
}
pkdensx = Vectorize( pkdensx1, "z" )
pkdensy1 = function( z, mean=0, sd=1 ) {
  h = bw.bcv( y )
  return( scdf1( z, y, h ) )
}
pkdensy = Vectorize( pkdensy1, "z" )

dkdensx1 = function( z, mean=0, sd=1 ) {
  h = bw.bcv( x )
  return( mean( dnorm( (z-x)/h )/h ) )
}
dkdensx = Vectorize( dkdensx1, "z" )
dkdensy1 = function( z, mean=0, sd=1 ) {
  h = bw.bcv( y )
  return( mean( dnorm( (z-y)/h )/h ) )
}
dkdensy = Vectorize( dkdensy1, "z" )

squantile = function( x, q, h = bw.ucv(x) ) {
  obj = function( t, x, h ) {
    mean( pnorm( (t-x)/h ) ) - q
  }
  objv = Vectorize( obj, "t" )
#  plot( seq(-5,5,0.01), objv( seq(-5,5,0.01), x=x,h=h ), type="l" )
  xmin = min( x )
  xmax = max( x )
  st = xmax - xmin
  while( obj(xmin,x,h) > 0 ) { xmin = xmin - st/10 }
  while( obj(xmax,x,h) < 0 ) { xmax = xmax + st/10 }
  uniroot( obj, c(xmin, xmax), tol = 0.0001, x = x, h = h )$root
}

squant = Vectorize( squantile, "q" )

qkdensx = function( q, mean=0, sd=1 ) {
  h = bw.bcv( x )
  return( squant( x, q, h ) )
}
qkdensy = function( q, mean=0, sd=1 ) {
  h = bw.bcv(y)
  return( squant( y, q, h ) )
}

# copula
library(copula)

normal.cop <- normalCopula( c(0), dim = 2 )
u <- apply( cbind(x,y), 2, rank) / (n + 1) ##
fit.tau <- fitCopula(normal.cop, u, method="itau")
fit.tau
# gives: param = 0.6338359 (used in mvdc call)

myMvd3 <- mvdc(copula = ellipCopula(family = "normal", param = 0.6338359),
margins = c("kdensx", "kdensy"), paramMargins = list(list(mean = mean(x),
sd = sd(x)), list(mean = mean(y), sd = sd(y))))

# now we can plot the copula

contour(myMvd3, dmvdc, xlim = c(-0.03, 0.03), ylim = c(-0.03, 0.03), las=1)

persp(myMvd3, dmvdc, xlim = c(-0.03, 0.03), ylim = c(-0.03, 0.03),  xlab = "DAX", ylab = "FTSE", zlab ="", las=1 )