R/Copulas.R

##########
# Defining Classes & Methods

##########
# Defining class copula

setClass("copula", representation(parameters = "numeric", param.names = "character", param.lowbnd = "numeric", param.upbnd = "numeric", message = "character"))

 ##########
 # Remove methods if they already exist

# if(isGeneric("pcopula")) removeGeneric("pcopula")
# if(isGeneric("dcopula")) removeGeneric("dcopula")
# if(isGeneric("rcopula")) removeGeneric("rcopula")
# if(isGeneric("contour.plot")) removeGeneric("contour.plot")

# if(isGeneric("CondCinverse")) removeGeneric("CondCinverse")

# if(isGeneric("A")) removeGeneric("A")
# if(isGeneric("AsecondDer")) removeGeneric("AsecondDer")
# if(isGeneric("AfirstDer")) removeGeneric("AfirstDer")

# if(isGeneric("Hderiv")) removeGeneric("Hderiv")
# if(isGeneric("PhiDer")) removeGeneric("PhiDer")
# if(isGeneric("Phi")) removeGeneric("Phi")
# if(isGeneric("InvPhi")) removeGeneric("InvPhi")
# if(isGeneric("InvPhisecondDer")) removeGeneric("InvPhisecondDer")
# if(isGeneric("InvPhifirstDer")) removeGeneric("InvPhifirstDer")
# if(isGeneric("DerPhi.minus1")) removeGeneric("DerPhi.minus1")
# if(isGeneric("tail.index")) removeGeneric("tail.index)
# if(isGeneric("dcdx")) removeGeneric("dcdx")
# if(isGeneric("lambda")) removeGeneric("lambda")

# if(isGeneric("Kendalls.tau")) removeGeneric("Kendalls.tau")
# if(isGeneric("Spearmans.rho")) removeGeneric("Spearmans.rho")

##############
# Defining class bivd

setClass("bivd", representation(Xmarg = "character", Ymarg = "character",   
     param.Xmarg = "numeric", param.Ymarg = "numeric"))


## if(isGeneric("pbivd")) removeGeneric("pbivd")
## if(isGeneric("dbivd")) removeGeneric("dbivd")
## if(isGeneric("rbivd")) removeGeneric("rbivd")


##############
# Defining class empirical copula

setClass("empirical.copula", representation(x = "numeric", y = "numeric"))


classNames <- c("normal","frank","kimeldorf.sampson",
                 "gumbel","galambos","husler.reiss","tawn",
                 "bb5","normal.mix",
                 "bb1","bb2","bb3","bb6","bb7","joe",
                 "bb4","empirical")


#############
# Set generic methods for class copula

setGeneric("pcopula", function(copula, u, v) standardGeneric("pcopula"))
setGeneric("dcopula", function(copula, u, v) standardGeneric("dcopula"))
setGeneric("rcopula", function(copula, n) standardGeneric("rcopula"))
setGeneric("contour.plot", function(object,n = 100, nlevels = 10, add = F,  ... ) standardGeneric("contour.plot"))

setGeneric("CondCinverse", function(copula, q, u ) standardGeneric("CondCinverse"))

setGeneric("A", function(copula, t) standardGeneric("A"))
setGeneric("AsecondDer", function(copula, t) standardGeneric("AsecondDer"))
setGeneric("AfirstDer", function(copula, t) standardGeneric("AfirstDer"))
setGeneric("Hderiv", function(copula, z) standardGeneric("Hderiv"))

setGeneric("Phi", function(copula, t) standardGeneric("Phi"))
setGeneric("PhiDer", function(copula, t) standardGeneric("PhiDer"))
setGeneric("InvPhi", function(copula, s) standardGeneric("InvPhi"))
setGeneric("InvPhisecondDer", function(copula, s) standardGeneric("InvPhisecondDer"))
setGeneric("InvPhifirstDer", function(copula, s) standardGeneric("InvPhifirstDer"))
setGeneric("DerPhi.minus1", function(copula, s) standardGeneric("DerPhi.minus1"))

setGeneric("Kendalls.tau", function(copula, tol = 1e-5) standardGeneric("Kendalls.tau"))
setGeneric("Spearmans.rho", function(copula, tol = 1e-5) standardGeneric("Spearmans.rho"))

setGeneric("lambda", function(copula, t) standardGeneric("lambda"))
setGeneric("tail.index", function(copula, ...) standardGeneric("tail.index"))
setGeneric("dcdx", function(copula, u, v) standardGeneric("dcdx"))


#############################################
## copula classes implemented
#############################################

 setClass("normal.copula","copula")
 setClass("frank.copula","copula")
 setClass("kimeldorf.sampson.copula","copula")

 setClass("gumbel.copula","copula")
 setClass("galambos.copula","copula")
 setClass("husler.reiss.copula","copula")
 setClass("tawn.copula","copula")
setClass("bb5.copula","copula")

setClass("normal.mix.copula","copula")

setClass("bb1.copula","copula")
setClass("bb2.copula","copula")
setClass("bb3.copula","copula")
setClass("bb4.copula","copula")
setClass("bb6.copula","copula")
setClass("bb7.copula","copula")
setClass("joe.copula","copula")


#############################################
## bivd classes implemented
#############################################

setClass("normal.bivd",representation("bivd", copula = "normal.copula"))
setClass("frank.bivd",representation("bivd", copula = "frank.copula"))
setClass("kimeldorf.sampson.bivd",representation("bivd", copula = "kimeldorf.sampson.copula"))

setClass("gumbel.bivd",representation("bivd", copula = "gumbel.copula"))
setClass("galambos.bivd",representation("bivd", copula = "galambos.copula"))
setClass("husler.reiss.bivd",representation("bivd", copula = "husler.reiss.copula"))
setClass("tawn.bivd",representation("bivd", copula = "tawn.copula"))
setClass("bb5.bivd",representation("bivd", copula = "bb5.copula"))
setClass("normal.mix.bivd",representation("bivd", copula = "normal.mix.copula"))

setClass("bb1.bivd",representation("bivd", copula = "bb1.copula"))
setClass("bb2.bivd",representation("bivd", copula = "bb2.copula"))
setClass("bb3.bivd",representation("bivd", copula = "bb3.copula"))
setClass("bb4.bivd",representation("bivd", copula = "bb4.copula"))
setClass("bb6.bivd",representation("bivd", copula =  "bb6.copula"))
setClass("bb7.bivd",representation("bivd", copula = "bb7.copula"))
setClass("joe.bivd",representation("bivd", copula = "joe.copula"))

setClass("empirical.bivd", representation(x = "numeric", y = "numeric"))

setMethod("show","copula", function(object) { cat("    ", object@message,".\n    Parameters : \n", sep = "")
     for (i in (1:length(object@parameters))) cat("      ", object@param.names[i]," = ",object@parameters[i],"\n")
   })


#############################################
## Prepare for method functions
#############################################

funNames <- c("pbivd","dbivd","rbivd","persp.dbivd","persp.pbivd",
               "contour.dbivd","contour.pbivd") 

funCall.pbivd <- "function(dist, x, y) "
nn.pbivd <- 4

funCall.dbivd <- "function(dist, x, y) "
nn.dbivd <- 4 

funCall.rbivd <-"function(dist,n) "
nn.rbivd <- 3

funCall.persp.dbivd <- "function(dist, n = 50, xlim = NA, ylim = NA, ...) "
funCall.persp.pbivd <- "function(dist, n = 50, xlim = NA, ylim = NA, ...) "
funCall.contour.dbivd <- "function(dist, n = 50, xlim = NA, ylim = NA, ...) "
funCall.contour.pbivd <- "function(dist, n = 50, xlim = NA, ylim = NA, ...) "

nn.persp.dbivd <- 6
nn.persp.pbivd <- 6
nn.contour.dbivd <- 6
nn.contour.pbivd <- 6


bivd <- function(cop, marginX = "unif", marginY = marginX, 
            param.marginX = c(0,1),
            param.marginY = param.marginX) {

     yyy <- class(cop)
       cl <- substring(yyy,1,nchar(yyy) - 7)

     val <- new( paste(cl,"bivd",sep = "."), copula = cop, Xmarg = marginX, Ymarg = marginY,
       param.Xmarg = param.marginX ,       
       param.Ymarg = param.marginY)

       val  
}

evalFunc <- function(x, fun, param) 
{
       n <- length(param)
       call <- switch(n,
                    expression(fun(x,param[1])),
                expression(fun(x,param[1],param[2])),
                 expression(fun(x,param[1],param[2],param[3])),
                 expression(fun(x,param[1],param[2],param[3],param[4])))
       eval(call)
}

pbivd <- function(dist, x, y) {
         xmar <- evalFunc(x,get(paste("p",dist@Xmarg,sep = "")),dist@param.Xmarg)
        ymar <- evalFunc(y,get(paste("p",dist@Ymarg,sep = "")),dist@param.Ymarg )
        return(pcopula(dist@copula,xmar,ymar))
}

dbivd <- function(dist, x, y) {
        xmar <- evalFunc(x,get(paste("p",dist@Xmarg,sep = "")),dist@param.Xmarg)
       ymar <- evalFunc(y,get(paste("p",dist@Ymarg,sep = "")),dist@param.Ymarg )
       dmarx <- evalFunc(x,get(paste("d",dist@Xmarg,sep = "")),dist@param.Xmarg)
       dmary <- evalFunc(y,get(paste("d",dist@Ymarg,sep = "")),dist@param.Ymarg)
       return(dcopula(dist@copula,xmar,ymar)*dmarx *dmary)
}

rbivd <- function(dist,n) {
   uv <- rcopula(dist@copula,n)
   Xsim <- evalFunc(uv$x, get(paste("q",dist@Xmarg,sep = "")),dist@param.Xmarg)
   Ysim <- evalFunc(uv$y, get(paste("q",dist@Ymarg,sep = "")),dist@param.Ymarg)
   val <- list(x = Xsim, y = Ysim)
   data.frame(val)
}

###########################################################################################################################################################################
########################################## 
### Functions pbivd
##########################################


#if(isGeneric("pbivd"))removeGeneric("pbivd")
setGeneric("pbivd", function(dist,...)standardGeneric("pbivd"))
setMethod("pbivd","ANY",function(dist, x, y)   {
    xmar <- evalFunc(x, get(paste("p", dist@Xmarg, sep = "")), dist@param.Xmarg)
    ymar <- evalFunc(y, get(paste("p", dist@Ymarg, sep = "")), dist@param.Ymarg)
    return(pcopula(dist@copula, xmar, ymar))
} )


## improved
setGeneric("dbivd",function(dist,...) standardGeneric("dbivd"))

setMethod("dbivd", "ANY",function(dist, x, y){
    xmar <- evalFunc(x, get(paste("p", dist@Xmarg, sep = "")), dist@param.Xmarg)
    ymar <- evalFunc(y, get(paste("p", dist@Ymarg, sep = "")), dist@param.Ymarg)
    dmarx <- evalFunc(x, get(paste("d", dist@Xmarg, sep = "")), dist@param.Xmarg)
    dmary <- evalFunc(y, get(paste("d", dist@Ymarg, sep = "")), dist@param.Ymarg)
    return(dcopula(dist@copula, xmar, ymar) * dmarx * dmary)
} )



# ###########################
### Functions rbivd  ### 
setGeneric("rbivd",function(dist,...) standardGeneric("rbivd"))
setMethod("rbivd","ANY",function(dist,n)  {
    uv <- rcopula(dist@copula, n)
    Xsim <- evalFunc(uv$x, get(paste("q", dist@Xmarg, sep = "")), dist@param.Xmarg)
    Ysim <- evalFunc(uv$y, get(paste("q", dist@Ymarg, sep = "")), dist@param.Ymarg)
    val <- list(x = Xsim, y = Ysim)
    data.frame(val)
} )


##############################################################################################################################
### Functions persp.dbivd  ### 

#if(isGeneric("persp.dbivd"))removeGeneric("persp.dbivd")
setGeneric("persp.dbivd",function(dist,...) standardGeneric("persp.dbivd"))
setMethod("persp.dbivd", "ANY",function(dist, n = 50, xlim = NA, ylim = NA, ...)   {
    divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
    divis <- divis[2:(n + 1)]
    if((is.na(xlim[1])) | (length(xlim) != 2)) {
        x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")), dist@param.Xmarg)
    }
    else {
        x <- seq(from = xlim[1], to = xlim[2], length = n)
    }
    if(is.na(ylim[1]) | length(ylim) != 2) {
        y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")), dist@param.Ymarg)
    }
    else {
        y <- seq(from = ylim[1], to = ylim[2], length = n)
    }
    xmat <- rep(x, n)
    ymat <- rep(y, each = n)
    zmat <- dbivd(dist, xmat, ymat)

persp3d(x, y, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "Density")

#    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat), z = as.vector(zmat)))
#    guiPlot("Data Grid Surface", DataSetValues = val2, ...)
    val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, data = zmat))
    invisible(val)
} )


 #
# 
### Functions persp.pbivd  ### 

setGeneric("persp.pbivd",function(dist,...) standardGeneric("persp.pbivd"))
setMethod("persp.pbivd","ANY",function(dist, n = 50, xlim = NA, ylim = NA, ...)  ## standardGeneric("persp.pbivd"))
## setMethod("persp.pbivd","normal.bivd",
##      function(dist, n = 50, xlim = NA, ylim = NA, ...)
          {
    divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
    divis <- divis[2:(n + 1)]
    if((is.na(xlim[1])) | (length(xlim) != 2)) {
        x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")), dist@param.Xmarg)
    }
    else {
        x <- seq(from = xlim[1], to = xlim[2], length = n)
    }
    if(is.na(ylim[1]) | length(ylim) != 2) {
        y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")), dist@param.Ymarg)
    }
    else {
        y <- seq(from = ylim[1], to = ylim[2], length = n)
    }
    xmat <- rep(x, n)
    ymat <- rep(y, each = n)
    zmat <- pbivd(dist, xmat, ymat)
#    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat), z = as.vector(zmat)))
#    guiPlot("Data Grid Surface", DataSetValues = val2, ...)

persp3d(x, y, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "Density")

     val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, data = zmat))
    invisible(val)
} )

### Functions contour.dbivd  ### 

#if(isGeneric("contour.dbivd")) removeGeneric("contour.dbivd")
setGeneric("contour.dbivd",function(dist,...) standardGeneric("contour.dbivd"))
setMethod("contour.dbivd","ANY",function(dist,n= 50, xlim = NA, ylim = NA, ...) #standardGeneric"contour.dbivd"))
## setMethod("contour.dbivd","normal.bivd",
##      function(dist, n = 50, xlim = NA, ylim = NA, ...) 
          {
    divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
    divis <- divis[2:(n + 1)]
    if((is.na(xlim[1])) | (length(xlim) != 2)) {
        x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")), dist@param.Xmarg)
    }
    else {
        x <- seq(from = xlim[1], to = xlim[2], length = n)
    }
    if(is.na(ylim[1]) | length(ylim) != 2) {
        y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")), dist@param.Ymarg)
    }
    else {
        y <- seq(from = ylim[1], to = ylim[2], length = n)
    }
    xmat <- rep(x, n)
    ymat <- rep(y, each = n)
    zmat <- dbivd(dist, xmat, ymat)
    val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, data = zmat))
contour(x,y,val$z)
title("Contour Plot of the Density", xlab="x",ylab="y")
    invisible(val)
} )



### Functions contour.pbivd  ### 

#if(isGeneric("contour.pbivd"))removeGeneric("contour.pbivd")
setGeneric("contour.pbivd",function(dist,...) standardGeneric("contour.pbivd"))
setMethod("contour.pbivd","ANY",function(dist, n = 50, xlim = NA, ylim = NA, ...) ##  standardGeneric("contour.pbivd"))


## setMethod("contour.pbivd","normal.bivd",
##      function(dist, n = 50, xlim = NA, ylim = NA, ...)
          {
    divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
    divis <- divis[2:(n + 1)]
    if((is.na(xlim[1])) | (length(xlim) != 2)) {
        x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")), dist@param.Xmarg)
    }
    else {
        x <- seq(from = xlim[1], to = xlim[2], length = n)
    }
    if(is.na(ylim[1]) | length(ylim) != 2) {
        y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")), dist@param.Ymarg)
    }
    else {
        y <- seq(from = ylim[1], to = ylim[2], length = n)
    }
    xmat <- rep(x, n)
    ymat <- rep(y, each = n)
    zmat <- pbivd(dist, xmat, ymat)
    val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, data = zmat))
contour(x,y,val$z)
title("Contour Plot of the CDF", xlab="x",ylab="y")
    invisible(val)
} )








##############################################################################################################################################################################

#########################
# Empirical  Copulas Methods


setMethod("Kendalls.tau","empirical.copula", function(copula, tol = 1e-5) 
{
   obs.x <- c(copula@x)
    obs.y <- c(copula@y)
    slistX <- sort.list(obs.x)
   slistY <- sort.list(obs.y)
   xsort <- as.double(obs.x[slistX])
   ysort <- as.double(obs.y[slistY])
   ind <- as.integer(match(slistX,slistY) - 1)
   nsamp <- as.integer(length(xsort))
   tau <- as.double(-999)

   outC <- .C("tau_EC", xsort,ysort,ind, nsamp,tau)
   tau <- outC[[5]]
   tau          
})  

setMethod("Spearmans.rho","empirical.copula", function(copula, tol = 1e-5) 
{
   obs.x <- c(copula@x)
    obs.y <- c(copula@y)
    slistX <- sort.list(obs.x)
   slistY <- sort.list(obs.y)
   xsort <- as.double(obs.x[slistX])
   ysort <- as.double(obs.y[slistY])
   ind <- as.integer(match(slistX,slistY) - 1)
   nsamp <- as.integer(length(xsort))
   tau <- as.double(-999)

   outC <- .C("rho_EC", xsort,ysort,ind, nsamp,tau)
   tau <- outC[[5]]
   tau          
}   )


setMethod("A","empirical.copula", function(copula, t) {

   u <- copula@x
   v <- copula@y

   p <- function (t) {
       1 - t
   }

    n <- length(u)
    N <- length(t)

    Z <- sort(log(u)/log(u*v))
    Q <- cumprod((Z/(1-Z))^(1/n))
    val <- vector(length = N, mode = "numeric")
    val[t < 0 | t > 1] <- 1

    block1 <- t >=0 & t <= Z[1]
    block2 <- t > Z[1] & t < Z[n]
    block3 <- t >= Z[n] & t <= 1


    Qn <- Q[n]
    val[block1] <- (1 - t[block1])*Qn^(1 - p(t[block1]))
    val[block3] <- (t[block3])*Qn^( - p(t[block3]))

    tt <- t[block2]
    nt <- length(tt)
    tsort <- sort(tt)
    N2 <- length(tsort)
    v2 <- vector(length = N2, mode = "numeric")
    v2 [1:N2] <- -1
    outC <- .C("empirfunc",as.double(tsort),as.double(Z),as.integer(N2),
           as.integer(n),as.integer(v2))[[5]]

    indback <- c(1: N2)
    indback[sort.list(tt)] <- c(1: N2)
    i <- outC[indback] 
    iovn <- i/n
    v2 <- (tt^(iovn))*((1-tt)^(1-iovn))*Qn^(1 - p(tt))/(Q[i])

    val[block2] <- v2

    #check fo consistency with 0.5 - 1 range. change to that range if need, but 
    #issue warning
    badval <- (val < 0.5 | val > 1)
    if (sum(badval) > 0) {
           warning(paste("Some of the estimated values of A(t) are not within [0.5,1] range",
           "They are chaged to fit this range", sep = "  "));
           val[val < 0.5] <- 0.5;
           val[val > 1] <- 1.0
        }

    val
})


################################################################################
# Fitting copula to data                                                       #
# Generic MLE method                                                           #
################################################################################

fit.copula <- function(data,family="normal", plotPicture=T, init.est=NA, epsilon=1e-5, ...)
{
     if (!inherits(data,"empirical.copula")) {data <- empirical.copula(data)}
     families.implemented <- c("normal","frank","kimeldorf.sampson",
                      "gumbel","galambos","husler.reiss",
                     "bb1","tawn","bb2","bb3","bb5","bb6","bb7","normal.mix","bb4", "joe")

     w <- charmatch(family, families.implemented, nomatch = -1)
     if (w <= 0) {
           message <- paste("Family ",family," is not implemented in fit.copula function")
           message <- paste(message," Valid family names are ",families.implemented)
           stop(message)
     }

     if (is.na(init.est)) {
     Kt <- Kendalls.tau(empirical.copula(data))
     if (Kt <= 0) {
           message <- paste("Negative dependence is detected in fit.copula function")
           message <- paste(message," Negative dependence is not implemented yet.")
           stop(message)
     }
    #starting point of parameters!
    Ktrd <- round(Kt, digits = 1)
    if (Ktrd == 0) {
          Ktrd <- 0.1
          message <- paste("Very small dependence is detected in fit.copula function")
          warning(message)
    }
    if (Ktrd == 1) Ktrd <- 0.9
    nline <- Ktrd *10 + 1
    x0 <- switch(w, Kend.table[nline,2],
          Kend.table[nline,3],Kend.table[nline,4],Kend.table[nline,5],Kend.table[nline,6],
          Kend.table[nline,7],  c(Kend.table[nline,9], Kend.table[nline,8]),
          c(Kend.table[nline,11],Kend.table[nline,11], Kend.table[nline,10]),
          c(Kend.table[nline,12],Kend.table[nline,13]),
          c(Kend.table[nline,14],Kend.table[nline,15]),
          c(Kend.table[nline,16],Kend.table[nline,17]),
          c(Kend.table[nline,18],Kend.table[nline,19]),
          c(Kend.table[nline,20],Kend.table[nline,21]),
          c(0.5, Kend.table[nline,2],Kend.table[nline,2]),
          c(Kend.table[nline,22],Kend.table[nline,23]),
          c(2.0) )
    }
    else{ x0 <- init.est }
    copula <- switch(w, normal.copula(x0), frank.copula(x0), kimeldorf.sampson.copula(x0),
            gumbel.copula(x0), galambos.copula(x0), husler.reiss.copula(x0),
            bb1.copula(x0[1],x0[2]), tawn.copula(x0[1],x0[2],x0[3]),
            bb2.copula(x0[1],x0[2]),
            bb3.copula(x0[1],x0[2]),
            bb5.copula(x0[1],x0[2]),
            bb6.copula(x0[1],x0[2]),
            bb7.copula(x0[1],x0[2]),
            normal.mix.copula(x0[1],x0[2],x0[3]),
            bb4.copula(x0[1],x0[2]), joe.copula(x0[1]))
    copTemp <- copula

    neg.log.lik <- function(param, copula, data) {
               copula@parameters <- param
               val <- dcopula(copula, data@x, data@y)
               nll <- sum(log(val))
               -nll
    }
    fit <- nlminb(x0, neg.log.lik, copula=copula, data=data, lower=copula@param.lowbnd, upper=copula@param.upbnd)
    cat("Maximum Likelihood Estimation of copula:\n")
    cat("MLE terminated due to ",fit$message,"\n")
    copula@parameters <- fit$par
    names(copula@parameters) <- names(copTemp@parameters)

    if (plotPicture) {
          contour.plot(data)
          contour.plot(copula, labex = 0, col = 5, lwd = 2, lty = 4, add = T, xlab = "", ylab = "")
          title("Empirical and Fitted Level Sets")  
    }
    copula
}



##########################
##########################
# Copula Families


######### Family B1 ######################
# Bivariate Normal                       #
##########################################

normal.copula <- function(delta) {
      if (is.na(delta) | delta > 1 | delta < 0)
          stop("Parameter of the bivariate the normal.copula should be between 0 and 1")

      val <- new("normal.copula", parameters = delta, param.names = "delta",
               param.lowbnd = 0, param.upbnd = 1, message = "Normal copula family")
     val
}


############################################################
# Mixtures of Normal copulas ###############################

normal.mix.copula <- function(p, delta1, delta2) {
      if (delta1 > 1 | delta1 < 0|
          p > 1 | p < 0 | delta2 > 1 | delta2 < 0)
          stop("All parameters for b1mix.copula should be between 0 and 1")

     val <- new("normal.mix.copula", parameters = c(p,delta1, delta2),
            param.names =  c("p","delta1", "delta2"),
             param.lowbnd = c(0,0,0), param.upbnd = c(1,1,1), 
             message = "Mixed Normal copula family")
    val    
}



###########################################################
#### Family B4 (Kimerdorf and Sampson 91975)) #############

kimeldorf.sampson.copula <- function(delta) {
   if (is.na(delta) | delta < 0)
          stop("Parameter delta for kimeldorf.sampson.copula should be >= 0")

   val <- new("kimeldorf.sampson.copula", parameters = delta,  param.names = "delta",     
               param.lowbnd = 0, param.upbnd = Inf, 
               message = "Kimeldorf and Sampson copula family") 
   val
}

setMethod("CondCinverse","kimeldorf.sampson.copula", function(copula, q, u) {
   delta <- as.vector(copula@parameters[1])
   ((q^(-delta/(delta +1)) -1)*u^(-delta)+1)^(-1/delta)
})

setMethod("Phi","kimeldorf.sampson.copula", function(copula, t) {
    delta <- as.vector(copula@parameters[1])
    t^(-delta) -1
})

setMethod("PhiDer","kimeldorf.sampson.copula", function(copula, t) {
     delta <- as.vector(copula@parameters[1])
     (-delta)* t^(-delta - 1)
})

##################################################################
##################   EXTREME Value copulas #######################
#
#NOTE: all extreme value copulas are 
#created in the archimax setting
# general method for generation for the EV copulas:


####################
# Family   B6  #####
####################

gumbel.copula <- function(delta) {
   if (is.na(delta) | delta < 1) {
      stop("For Gumbel family parameter delta >= 1")
   }
   val <- new("gumbel.copula", parameters = delta, param.names = "delta",
              param.lowbnd = 1, param.upbnd = Inf,
              message = "Gumbel copula family; Extreme value copula")
    val
}

setMethod("A","gumbel.copula", function(copula, t) {
  delta <- as.vector(copula@parameters[1])
  (t^delta+(1-t)^delta)^(1/delta)
})

setMethod("Hderiv","gumbel.copula", function(copula, z) {
   d <- as.vector(copula@parameters[1])
  (d*(-((-1 + z)*z))^(-1 + d))/((1 - z)^d + z^d)^2
})

setMethod("AfirstDer","gumbel.copula", function(copula, t) {
    delta <- as.vector(copula@parameters[1])
        ((-(delta*(1 - t)^(-1 + delta)) + delta*t^(-1 + delta))*
    ((1 - t)^delta + t^delta)^(-1 + delta^(-1)))/delta
    })

setMethod("AsecondDer","gumbel.copula",function(copula, t) {
   d <-  as.vector(copula@parameters[1])
   (-1 + d)*(-((-1 + t)*t))^(-2 + d)*((1 - t)^d + t^d)^
      (-2 + d^(-1))
})

######################################
## Family b7                         #
######################################

galambos.copula <- function(delta) {
   if (is.na(delta) | delta < 0) {
      stop("Parameter delta >=  0 for B7 copula family") }

   val <- new("galambos.copula", parameters = delta, param.names = "delta", 
         param.lowbnd = 0, param.upbnd = Inf,
         message = "Galambos copula family; Extreme value copula")

   val  
}

setMethod("A","galambos.copula", function(copula, t) { 
  delta <- as.vector(copula@parameters[1])
  1 - (t^(-delta)+(1-t)^(-delta))^(-1/delta)
})

setMethod("AfirstDer","galambos.copula", function(copula,t) {
   delta <- as.vector(copula@parameters[1])
  ((1 - t)^(-1 - delta) - t^(-1 - delta))/
    ((1 - t)^(-delta) + t^(-delta))^((1 + delta)/delta)
})

setMethod("AsecondDer","galambos.copula", function(copula, t) {
   delta <-  as.vector(copula@parameters[1])
((1 + delta)*(-((-1 + t)*t))^(-2 + delta))/
  (((1 - t)^(-delta) + t^(-delta))^delta^(-1)*
    ((1 - t)^delta + t^delta)^2)
})

######################################
## Family b8                         #
######################################

husler.reiss.copula <- function(delta) {

    if (is.na(delta) | delta < 0) {
      stop("Parameter delta >=  0 for Husler and Reiss copula family") }

  val <- new("husler.reiss.copula", parameters = delta, param.names = "delta",
             param.lowbnd = 0, param.upbnd = Inf,
           message = "Husler and Reiss copula family; Extreme value copula")
    val 

}

setMethod("A","husler.reiss.copula", function(copula, t) {

  d <- as.vector(copula@parameters[1])  
  ft1 <- pnorm(1/d + 1/2*d*log(t/(1 - t)))
  ft2 <- pnorm(1/d + 1/2*d*log((1-t)/t))

  t*ft1 + (1 - t) *ft2
})

setMethod("AfirstDer","husler.reiss.copula", function(copula, t) {
    d <- as.vector(copula@parameters[1])
   ft1 <- pnorm(1/d + 1/2*d*log(t/(1 - t)))
   ft2 <- pnorm(1/d + 1/2*d*log((1-t)/t))

   dn1 <- dnorm(1/d + 1/2*d*log(t/(1 - t))) * d/2/(1 - t)
   dn2 <- dnorm(1/d - 1/2*d*log(t/(1 - t)))  *d/2/t

    dn1 + ft1 - dn2 - ft2

})

setMethod("AsecondDer","husler.reiss.copula",function(copula, t) {
   d <- as.vector(copula@parameters[1])
   x1 <- 1/d + 1/2*d*log(t/(1 - t))
   x2 <- 1/d - 1/2*d*log(t/(1 - t))
   dterm <- d/2/t/(1-t)

  dnorm(x1)*dterm/(1-t)*(1 - x1*d/2) +
    dnorm(x2)*dterm/t*(1 - x2*d/2)
})


######################################
## Family Tawn                       #
######################################

tawn.copula <- function(a,b,r) 
{
   if (a < 0 | r < 1 |b > 1)
      stop("Parameter a and b >=  0 for Tawn family")

    val <- new("tawn.copula", parameters = c(a, b, r), param.names = c("alpha", "lambda", "r"),
              param.lowbnd = c(0,0,1), param.upbnd = c(1, 1, Inf),  
            message = "Tawn copula family; Extreme Value copula")
    val  
}


setMethod("A","tawn.copula", function(copula, t) {
  a <- as.vector(copula@parameters[1])
  b <- as.vector(copula@parameters[2])
  r <- as.vector(copula@parameters[3])
  1 - b + (b - a)*t + ((a* t)^r + (b*(1 - t))^r)^(1/r)
})

setMethod("AfirstDer","tawn.copula", function(copula, t) 
{
  a <- as.vector(copula@parameters[1])
  b <- as.vector(copula@parameters[2])
  r <- as.vector(copula@parameters[3])
-a + b + ((-(b*r*(b*(1 - t))^(-1 + r)) +
       a*r*(a*t)^(-1 + r))*
     ((b*(1 - t))^r + (a*t)^r)^(-1 + r^(-1)))/r
})

setMethod("Hderiv","tawn.copula",function(copula, z) 
{
   a <- as.vector(copula@parameters[1])
  b <- as.vector(copula@parameters[2])
  r <- as.vector(copula@parameters[3])
1 + ((1 - 2*z)*(-a + b +
    ((-(b^r*r*(1 - z)^(-1 + r)) + a^r*r*z^(-1 + r))*
      (b^r*(1 - z)^r + a^r*z^r)^(-1 + r^(-1)))/r))/
  (1 + b*(-1 + z) - a*z + (b^r*(1 - z)^r + a^r*z^r)^
    r^(-1)) +
 ((1 - z)*z*
   (-(-a + b + ((-(b^r*r*(1 - z)^(-1 + r)) +
          a^r*r*z^(-1 + r))*(b^r*(1 - z)^r + a^r*z^r)^
          (-1 + r^(-1)))/r)^2 + a^r*b^r*(-1 + r)*
     (-((-1 + z)*z))^(-2 + r)*(b^r*(1 - z)^r + a^r*z^r)^
      (-2 + r^(-1))*(1 + b*(-1 + z) - a*z +
      (b^r*(1 - z)^r + a^r*z^r)^r^(-1))))/
  (1 + b*(-1 + z) - a*z + (b^r*(1 - z)^r + a^r*z^r)^r^(-1))^2
})

setMethod("AsecondDer","tawn.copula", function(copula, t) 
{ 
  a <- as.vector(copula@parameters[1])
  b <- as.vector(copula@parameters[2])
  r <- as.vector(copula@parameters[3])
  a^r*b^r*(-1 + r)*(-((-1 + t)*t))^(-2 + r)*
     (b^r*(1 - t)^r + a^r*t^r)^(-2 + r^(-1))
})


###########################################################
#### Family B3 (Frank (1979)) #############################

frank.copula <- function(delta) {
    if (is.na(delta) | delta < 0)
          stop("Parameter delta for frank.copula should be >= 0")

   val <- new("frank.copula", parameters = delta,  param.names = "delta",     
               param.lowbnd = 0, param.upbnd = Inf, 
               message = "Frank copula family") 
   val
}

setMethod("CondCinverse","frank.copula", function(copula, q, u) {
   delta <- as.vector(copula@parameters[1])
   eta <-  1 - exp(-delta)
   -1/delta*log(1-eta/((1/q - 1)*exp(-delta*u) + 1))
})

setMethod("Phi","frank.copula", function(copula, t) {
    delta <- as.vector(copula@parameters[1])
    -log((exp(-delta*t) - 1)/(exp(-delta) -1))
})


setMethod("PhiDer","frank.copula", function(copula, t) {
     delta <- as.vector(copula@parameters[1])
     delta/(exp((delta*t))*(-1 + exp(-(delta*t))))
})

############################################
#############BB1 copula ####################

archm.copula <- function(family = "BB1", param = c(0.8,1.43)) {
     family <- casefold(family, upper=F)

     families.implemented <- c("gumbel","frank",
         "kimeldorf.sampson","joe","bb1","bb2","bb3","bb6","bb7")

       w <- charmatch(family, families.implemented, nomatch = -1)

       if (w <= 0) {
           message <- paste("Family ",family," is not implemented in function")
        message <- paste(message," Valid family names are ",families.implemented)
            print(message)
           stop(message)
        }

       copula <- switch(w, gumbel.copula(param), frank.copula(param), 
               joe.copula(param), bb1.copula(param[1],param[2]),
               bb2.copula(param[1],param[2]),
               bb3.copula(param[1],param[2]),
              bb6.copula(param[1],param[2]),
              bb7.copula(param[1],param[2]))
       copula   

}

bb1.copula <- function(theta,delta) {
    if (is.na(delta) | is.na(theta) | delta < 1 | theta < 0)
      stop("Parameter delta >= 1 and theta >  0 for BB1 family")
   val <- new("bb1.copula", parameters = c(theta,delta),
            param.names = c("theta","delta"),
            param.lowbnd = c(0,1), param.upbnd = c(Inf, Inf),
   message = "BB1 copula family; Archimedean Copula")
   val
}

setMethod("Phi","bb1.copula", function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   delta <- as.vector(copula@parameters[2])
   (t ^(-theta) - 1)^(delta)
})


setMethod("PhiDer","bb1.copula",function(copula,t) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
     - delta*theta*(t ^(-theta) - 1)^(delta - 1)* t ^(-theta-1)
})

setMethod("InvPhi","bb1.copula",function(copula,s) {
     theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
        ( 1 + s^(1/delta) ) ^(-1/theta)
})

setMethod("InvPhifirstDer","bb1.copula", function(copula,s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
        (-1/theta)*( 1 + s^(1/delta) ) ^(-1/theta -1)*(1/delta)*s^(1/delta-1)
})

setMethod("InvPhisecondDer","bb1.copula", function(copula, s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
      (s^(-2 + delta^(-1))*(1 + s^delta^(-1))^
     (-2 - theta^(-1))*
    ((-1 + delta)*theta +
      s^delta^(-1)*(1 + delta*theta)))/(delta^2*theta^2)
})


###################################################################
###  Family   BB2                                               ###
###################################################################


bb2.copula <- function(theta,delta) {
    if (is.na(delta) | is.na(theta) | delta < 0 | theta < 0)
      stop("Parameter delta >= 0 and theta >  0 for BB2 family")

    val <- new("bb2.copula", parameters = c(theta,delta), 
            param.names = c("theta","delta"),
            param.lowbnd = c(0,0), param.upbnd = c(Inf, Inf), 
            message = "BB2 copula family; Archimedean Copula") 
      val
}

setMethod("Phi","bb2.copula", function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   delta <- as.vector(copula@parameters[2])
   -1 + exp(delta*(-1 + t^(-theta)))
})

setMethod("PhiDer","bb2.copula", function(copula,t) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
    -(delta*exp(delta*(-1 + t^(-theta)))*t^(-1 - theta)*theta)
})

setMethod("InvPhi","bb2.copula", function(copula,s) {
     theta <- as.vector(copula@parameters[1])
     delta <- as.vector(copula@parameters[2])
  (1 + log(1 + s)/delta)^(-theta^(-1))
})

setMethod("InvPhisecondDer","bb2.copula", function(copula, s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
     (1 + theta + delta*theta + theta*log(1 + s))/
  ((1 + s)^2*theta^2*(delta + log(1 + s))^2*
    ((delta + log(1 + s))/delta)^theta^(-1))
})


####################################################
## Family BB3  #####################################
####################################################

bb3.copula <- function(theta,delta) {
    if (is.na(delta) | is.na(theta) | delta < 0 | theta < 1)
      stop("Parameter delta >= 0 and theta >  1 for BB3 family")

    val <- new("bb3.copula", parameters = c(theta,delta), 
            param.names = c("theta","delta"),
       param.lowbnd = c(1,0), param.upbnd = c(Inf, Inf), 
     message = "BB3 copula family; Archimedean Copula") 
    val
}

setMethod("Phi","bb3.copula",function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   delta <- as.vector(copula@parameters[2])
    -1 + exp(delta*(-log(t))^theta)

 })

setMethod("PhiDer","bb3.copula",function(copula,t) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
    (delta*exp(delta*(-log(t))^theta)*theta*(-log(t))^theta)/
                   (t*log(t))
})

setMethod("InvPhi","bb3.copula", function(copula,s) {
     theta <- as.vector(copula@parameters[1])
     delta <- as.vector(copula@parameters[2])
   exp(-(log(1 + s)/delta)^theta^(-1))
})

setMethod("InvPhisecondDer","bb3.copula",  function(copula, s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
  ((log(1 + s)/delta)^theta^(-1)*
    (-1 + theta + theta*log(1 + s) +
      (log(1 + s)/delta)^theta^(-1)))/
  (exp(log(1 + s)/delta)^theta^(-1)*(1 + s)^2*theta^2*
    log(1 + s)^2)
})


###################################################################
# Archimax Copulas #################################################
###################################################################

setMethod("Hderiv","bb4.copula", function(copula, z) {
       ad0 <- A(copula, z)
       ad1 <- AfirstDer(copula, z)
       ad2 <- AsecondDer(copula, z)

       1 + (1 - 2*z)*ad1/ad0 + z*(1 - z) *(ad2/ad0 - ad1^2/ad0^2)
    }   )

#################################################################
###   BB4 copula                                              ###
################################################################# 

bb4.copula <- function(theta,delta) {
    if (is.na(delta) | is.na(theta) | delta < 0 | theta < 0)
      stop("Parameter delta > 0 and theta >=  0 for BB4 family")

   val <- new("bb4.copula", parameters = c(theta,delta), 
            param.names = c("theta","delta"),
     param.lowbnd = c(0,0), param.upbnd = c(Inf, Inf), 
     message = "BB4 copula family; Archimax copula") 

   val
}

setMethod("Phi","bb4.copula", function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   t^(-theta) - 1
})

setMethod("PhiDer","bb4.copula",function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   (-theta)*t^(-theta - 1)
})

setMethod("InvPhi","bb4.copula",function(copula,s) {
   theta <- as.vector(copula@parameters[1])
   (s + 1)^(-1/theta)
})

setMethod("InvPhifirstDer","bb4.copula",function(copula,s) {
   theta <- as.vector(copula@parameters[1])
   (-1/theta)*(s + 1)^(-1/theta - 1)
})

setMethod("InvPhisecondDer","bb4.copula",function(copula,s) {
   theta <- as.vector(copula@parameters[1])
   (-1/theta)*(-1/theta - 1)*(s + 1)^(-1/theta - 2)
})

setMethod("DerPhi.minus1","bb4.copula",function(copula,s) {
   theta <- as.vector(copula@parameters[1])
    (-(s/theta))^(-1 - theta)^(-1)
})

setMethod("A","bb4.copula", function(copula, t) {
    delta <- as.vector(copula@parameters[2])
    1 - (t^(-delta)+(1-t)^(-delta))^(-1/delta)
})

setMethod("AfirstDer","bb4.copula", function(copula,t) {
   delta <- as.vector(copula@parameters[2])
   ((1 - t)^(-1 - delta) - t^(-1 - delta))/
    ((1 - t)^(-delta) + t^(-delta))^((1 + delta)/delta)
})

setMethod("AsecondDer","bb4.copula", function(copula, t) {
   delta <-  as.vector(copula@parameters[2])
 ((1 + delta)*(-((-1 + t)*t))^(-2 + delta))/
  (((1 - t)^(-delta) + t^(-delta))^delta^(-1)*
    ((1 - t)^delta + t^delta)^2)
})

###################################################################
### Family BB5  ##################################################


bb5.copula <- function(theta,delta) {
   if (is.na(delta) | delta < 0 | theta < 1) {
      stop("Parameter delta >  0 and theta >= 1 for BB5 copula family") }

   val <- new("bb5.copula", parameters = c(theta,delta), 
            param.names = c("theta","delta"),
            param.lowbnd = c(1,0), param.upbnd = c(Inf,Inf),
         message = "BB5 copula family; Extreme value copula")
    val  
}


setMethod("A","bb5.copula", function(copula, t) { 
  theta <- as.vector(copula@parameters[1])
  delta <- as.vector(copula@parameters[2])

  ((1 - t)^theta + t^theta - 
    ((1 - t)^(-(delta*theta)) + t^(-(delta*theta)))^
     (-delta^(-1)))^theta^(-1)
})

setMethod("AfirstDer","bb5.copula" ,function(copula,t) {
  theta <- as.vector(copula@parameters[1])
  delta <- as.vector(copula@parameters[2])
  ((1 - t)^theta + t^theta -
     ((1 - t)^(-(delta*theta)) + t^(-(delta*theta)))^
      (-delta^(-1)))^(-1 + theta^(-1))*
  (-(1 - t)^(-1 + theta) + t^(-1 + theta) +
    ((1 - t)^(-1 - delta*theta) - t^(-1 - delta*theta))/
     ((1 - t)^(-(delta*theta)) + t^(-(delta*theta)))^
      ((1 + delta)/delta))

})


setMethod("AsecondDer","bb5.copula", function(copula, t) {
   theta <- as.vector(copula@parameters[1])
  delta <- as.vector(copula@parameters[2])

((1 - t)^theta + t^theta - 
     ((1 - t)^(-(delta*theta)) + t^(-(delta*theta)))^
      (-delta^(-1)))^(-2 + theta^(-1))*
  ((-(1 - t)^(-1 + theta) + t^(-1 + theta) + 
        ((1 - t)^(-1 - delta*theta) - 
           t^(-1 - delta*theta))/
         ((1 - t)^(-(delta*theta)) + 
            t^(-(delta*theta)))^((1 + delta)/delta))^2*
     (-1 + theta^(-1))*theta + 
    ((1 - t)^theta + t^theta - 
       ((1 - t)^(-(delta*theta)) + t^(-(delta*theta)))^
        (-delta^(-1)))*((1 - t)^(-2 + theta)*
        (-1 + theta) + t^(-2 + theta)*(-1 + theta) - 
       ((1 + delta)*(-(1 - t)^(delta*theta) + 
             (1 - t)^(delta*theta)*t + 
             t^(1 + delta*theta))^2*theta)/
        ((-1 + t)^2*t^2*
          ((1 - t)^(-(delta*theta)) + 
             t^(-(delta*theta)))^delta^(-1)*
          ((1 - t)^(delta*theta) + t^(delta*theta))^2) + 
       ((-(1 - t)^(-2 - delta*theta) - 
            t^(-2 - delta*theta))*(-1 - delta*theta))/
        ((1 - t)^(-(delta*theta)) + t^(-(delta*theta)))^
         ((1 + delta)/delta)))

})


####################################################
## Family BB6  #####################################

bb6.copula <- function(theta,delta) {
    if (is.na(delta) | is.na(theta) | delta < 1 | theta < 1)
      stop("Parameter delta >= 1 and theta >=  1 for BB6 family")

    val <- new("bb6.copula", parameters = c(theta,delta), 
            param.names = c("theta","delta"),
            param.lowbnd = c(1,1), param.upbnd = c(Inf, Inf), 
            message = "BB6 copula family; Archimedean Copula") 
     val
}

setMethod("Phi","bb6.copula", function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   delta <- as.vector(copula@parameters[2])
 (-log(1 - (1 - t)^theta))^delta
 })

setMethod("PhiDer","bb6.copula",function(copula,t) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
    -((delta*(1 - t)^(-1 + theta)*theta*
      (-log(1 - (1 - t)^theta))^(-1 + delta))/
    (1 - (1 - t)^theta))
})

setMethod("InvPhi","bb6.copula", function(copula,s) {
     theta <- as.vector(copula@parameters[1])
     delta <- as.vector(copula@parameters[2])
   1 - (1 - exp(-s^(1/delta)))^(1/theta)
})

setMethod("InvPhisecondDer","bb6.copula",function(copula, s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])

    s1d <- exp(-s^delta^(-1))

  ((1 - s1d)^theta^(-1)*
    s^(-2 + delta^(-1))*
    ((-1 + delta)*(-1 + s1d)*theta + 
      s^delta^(-1)*(-1 + theta/s1d)))/
  (delta^2*(-1 + 1/s1d)^2*theta^2)
})

#####################################################################
####### Family BB7   ################################################

bb7.copula <- function(theta,delta) {
    if (is.na(delta) | is.na(theta) | delta < 0 | theta < 1)
      stop("Parameter delta > 0 and theta >=  1 for BB7 family")

   val <- new("bb7.copula", parameters = c(theta,delta), 
            param.names = c("theta","delta"),
            param.lowbnd = c(1,0), param.upbnd = c(Inf, Inf), 
            message = "BB7 copula family; Archimedean Copula") 
   val
}

setMethod("Phi","bb7.copula", function(copula,t) {
   theta <- as.vector(copula@parameters[1])
   delta <- as.vector(copula@parameters[2])
   -1 + (1 - (1 - t)^theta)^(-delta)
})

setMethod("PhiDer","bb7.copula", function(copula,t) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
   -(delta*(1 - (1 - t)^theta)^(-1 - delta)*
    (1 - t)^(-1 + theta)*theta)
})

setMethod("InvPhi","bb7.copula", function(copula,s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
        1 - (1 - (1 + s)^(- 1/delta))^( 1/theta)
})

setMethod("InvPhifirstDer","bb7.copula", function(copula,s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
      - ( 1/theta)* (1 - (1 + s)^(- 1/delta))^( 1/theta - 1)*(1 + s)^(- 1/delta - 1)*(1/delta)
})

setMethod("InvPhisecondDer","bb7.copula", function(copula, s) {
     theta <- as.vector(copula@parameters[1])
        delta <- as.vector(copula@parameters[2])
((1 - (1 + s)^(-delta^(-1)))^theta^(-1)*
    (-1 + (1 + s)^delta^(-1)*theta +
      delta*(-1 + (1 + s)^delta^(-1))*theta))/
  (delta^2*(1 + s)^2*(-1 + (1 + s)^delta^(-1))^2*
    theta^2)
})

###################################################################
###   Joe Copula ##################################################

joe.copula <- function(theta) {
    if (  is.na(theta) | theta < 1)
      stop("Parameter delta >= 1 for Joe family")

   val <- new("joe.copula", parameters = c(theta), param.names = "theta",
       param.lowbnd = c(1), param.upbnd = c(Inf), 
         message = "Joe copula family; Archimedean Copula") 
   val
}

setMethod("Phi","joe.copula", function(copula,t) {
   theta <- as.vector(copula@parameters[1])
    (-log(1 - (1 - t)^theta))
 })

setMethod("PhiDer","joe.copula", function(copula,t) {
         theta <- as.vector(copula@parameters[1])
            -((1 - t)^(-1 + theta)*theta)/(1 - (1 - t)^theta)
})

setMethod("InvPhi","joe.copula", function(copula,s) {
    theta <- as.vector(copula@parameters[1])
    1 - (1 - exp(-s))^(1/theta)
})

setMethod("InvPhisecondDer","joe.copula", function(copula, s) {
     theta <- as.vector(copula@parameters[1])
     s1d <- exp(-s)

  ((1 - s1d)^theta^(-1)*
    s^(-2 + 1)*
    ( s*(-1 + theta/s1d)))/((-1 + 1/s1d)^2*theta^2)
})

###################
# Tail Index Computations

setMethod("tail.index","bb1.copula", function(copula,...) 
{   
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
   LI <- 2^(-1/delta/theta)
    UI <- 2 - 2^(1/delta)
   val <- c(LI,UI)
   names(val) <- c("lower.tail","upper.tail")
   val
})

setMethod("tail.index","bb2.copula", function(copula,...) 
{
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
   LI <- 1
    UI <- 0
   val <- c(LI,UI)
   names(val) <- c("lower.tail","upper.tail")
   val
})

 setMethod("tail.index","bb3.copula", function(copula,...) 
{   
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
   if(theta == 1) {LI <- 2^(-1/delta)}
   else{LI <- 1}
    UI <- 2 - 2^(1/theta)
   val <- c(LI,UI)
   names(val) <- c("lower.tail","upper.tail")
   val
})

setMethod("tail.index","bb4.copula", function(copula,...) 
{   
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
   LI <- (2 - 2^(-1/delta))^(-1/theta)
   UI <- 2^(-1/delta)
   val <- c(LI,UI)
   names(val) <- c("lower.tail","upper.tail")
   val
})

setMethod("tail.index","bb6.copula", function(copula,...) 
{   
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
   LI <- 0
   UI <- 2 - 2^(1/delta/theta)
   val <- c(LI,UI)
   names(val) <- c("lower.tail","upper.tail")
   val
})

setMethod("tail.index","bb7.copula", function(copula,...) 
{   
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
   LI <- 2^(-1/delta)
   UI <- 2 - 2^(1/theta)
   val <- c(LI,UI)
   names(val) <- c("lower.tail","upper.tail")
   val
  })

setMethod("tail.index","empirical.copula", function(copula,...) 
{   
      n <- length(cop$x)
      alpha <- c(trunc(n/5):n)/n
      uind <- (1 - 2*alpha + pcopula(cop,alpha,alpha))/(1 - alpha)
      lind <- (pcopula(cop,1 - alpha,1 - alpha))/(1 - alpha)

      par(mfrow = c(2,1))
      plot(alpha,lind, ylab = "Lower Index")
      plot(1 - alpha,uind, ylab = "Upper Index", xlab = "alpha")
      par(mfrow = c(1,1))     
})


#################
# Functions for Comparing Multiple Copulas

cop.dist.L2.norm <- function(copula1, copula2, n = 100) {

    xmat <- rep(c(1:(n-1))/(n), n -1)
    ymat <- rep(c(1:(n-1))/n, each = n - 1) 
    zmat1 <- matrix(ncol = n - 1, nrow = n - 1, pcopula(copula1, xmat, ymat), byrow = F)    
    zmat2 <- matrix(ncol = n - 1, nrow = n - 1, pcopula(copula2, xmat, ymat), byrow = F)
    val <- sum((zmat2 - zmat1)^2)/n/n
    val 
}

A.dist.L1.norm <- function(copula1, copula2, n = 100) {

    inputPOS1 <- inherits(copula1,"empirical.copula") & 
                          hasMethod("A",class(copula2))
    inputPOS2 <- inherits(copula2,"empirical.copula") & 
                         hasMethod("A",class(copula1))
    if (!(inputPOS1 |  inputPOS2)) {
          warning("One of the copulas should be empirical copula, and another should be EV copula")
          return(NA)
    }

    t <- c(0:100)/100

    a1 <- A(copula1,t)
    a2 <- A(copula2,t)
    val <- sum(abs(a1 - a2))/100
    val 
}

Pearson.test <- function(empir.cop, param.copula, Ncells = 36, est.param = 0) 
{
    n <- floor(sqrt(Ncells + 1))
    data <- list(x = empir.cop@x, y = empir.cop@y)
    xll <- c(0:(n-1))/n
    xur <- c(1:n)/n
    yll <- c(0:(n-1))/n
    yur <- c(1:n)/n

    obs.x <- data$x
    obs.y <- data$y
    slistX <- sort.list(obs.x)
    slistY <- sort.list(obs.y)
    xsort <- as.double(obs.x[slistX])
    ysort <- as.double(obs.y[slistX])

    ind <- as.integer(match(slistX, slistY) - 1)
   N <- as.integer(n)
    nsamp <- as.integer(length(xsort))
        val1 <- as.double(rep(0, N*N))
    outC <- .C("jointDensity",
        xsort,
        ysort,
        nsamp,
       as.double(xur), as.double(yur),
        N,N,val1,as.double(1/n),as.double(1/n))
    val1 <- outC[[8]]
    emp.mat <- matrix(nrow = n, ncol = n, data = val1, byrow = T)

    copula <- param.copula
    x <- rep(xll,n)
    y <- rep(yll,each = n)
    v00 <- pcopula(copula,x,y)
    v00[ x==0 | y == 0] <- 0
    v00[ x == 1] <- y[x ==1]
    v00[ y == 1] <- x[y ==1]

    x <- rep(xur,n)
    y <- rep(yll,each = n)
    v01 <- pcopula(copula,x,y)
    v01[ x==0 | y == 0] <- 0
    v01[ x == 1] <- y[x ==1]
    v01[ y == 1] <- x[y ==1]

    x <- rep(xll,n)
    y <- rep(yur,each = n)
    v10 <- pcopula(copula,x,y)
    v10[ x==0 | y == 0] <- 0
    v10[ x == 1] <- y[x ==1]
    v10[ y == 1] <- x[y ==1]

    x <- rep(xur,n)
    y <- rep(yur,each = n)
    v11 <- pcopula(copula,x,y)
    v11[ x==0 | y == 0] <- 0
    v11[ x == 1] <- y[x ==1]
    v11[ y == 1] <- x[y ==1]

    vcop <-  matrix(nrow = n, ncol = n, data = (v11 - v01 -v10 +v00), byrow = T)

    #pull together all cells with p < 0.05!
    small.cells <- vcop < 0.005
    total.pr <- sum(vcop[small.cells])
    total.count <- sum(emp.mat[small.cells])

    #cat(sum(small.cells), total.pr)
    vcop[small.cells] <- 0
    emp.mat[small.cells] <- 0

    real.mat <-  vcop*nsamp     

    Chi <- sum((real.mat[!small.cells] - emp.mat[!small.cells])^2/real.mat[!small.cells])          

    if (total.pr > 0) Chi <- Chi + (total.count - total.pr *nsamp)^2/total.pr/nsamp

    Df <- as.integer(n*n) - 1 - est.param

    res <- list(Chi.Statistics = Chi, deegrees.freedom = Df, p.value = 1 - pchisq(Chi,Df))
    as.data.frame(res)
}

############################
# K estimators for copulas
##KEST

Wis <- function(x,y) {

   tf <- function(i,x,y) {
      tt <- x < x[i] & y < y[i]
      sum(tt)
   }

   ww <- matrix(ncol = length(x), nrow = 1, data = c(1:length(x)))
   ww <- apply(ww,2,tf,x = x, y = y)

  as.vector(ww)/(length(x) -1)
}

#lambda.archm.copula <- function(copula, t) {
#       Phi(copula,t)/PhiDer(copula,t)
#   }


K.Kern.est <- function(w,t, D) {
       b <- 1/length(w)*D
       tf <- function(t,w,b) {
           sum(pnorm(t-w,0, b))
       }

         ww <- matrix(ncol = length(t), nrow = 1, data = t)
        ww <- apply(ww,2,tf,w = w, b = b)
        ww[t ==0] <- 0
        ww[t ==1] <- length(w)         
      as.vector(ww)/(length(w))
    }

K.ast.Kern.est <- function(w,t, D) {
       b <- 1/length(w)*D
       tf <- function(t,w,b) {
           sum(pnorm(t-w,0, b))
       }

         ww <- matrix(ncol = length(t), nrow = 1, data = t)
        ww <- apply(ww,2,tf,w = w, b = b)
        ww[t ==0] <- 0
        ww[t ==0] <- 1
        val <- as.vector(ww)/(length(w))
        cummax(val)
    }

setMethod("lambda","empirical.copula", function(copula, t) {
         w <- Wis(copula@x, copula@y)
         t -  K.Kern.est(w,t, 0.5)
})




########################## 
### Functions dcopula  ### 


setMethod("dcopula","normal.copula",  function(copula, u, v) {
       delta <- copula@parameters[1]
      x <- qnorm(u)
       y <- qnorm(v)
       val <- (1 - delta^2)^(-1/2)*exp(-1/2*(1 - delta^2)^(-1)*
                      (x^2 + y^2 - 2* delta*x*y))*exp(1/2*(x^2 + y^2))
       val[ u <= 0 | v <= 0 ] <- 0
       val[ u >= 1 | v >= 1] <- 0
      val
} )


setMethod("dcopula","normal.mix.copula", function(copula, u, v) {
    p <- as.vector(copula@parameters[1])
    delta1 <- as.vector(copula@parameters[2])
   delta2 <- as.vector(copula@parameters[3])

      x <- qnorm(u)
       y <- qnorm(v)
       val1 <- (1 - delta1^2)^(-1/2)*exp(-1/2*(1 - delta1^2)^(-1)*
                      (x^2 + y^2 - 2* delta1*x*y))*exp(1/2*(x^2 + y^2)) 

       val2 <- (1 - delta2^2)^(-1/2)*exp(-1/2*(1 - delta2^2)^(-1)*
                      (x^2 + y^2 - 2* delta2*x*y))*exp(1/2*(x^2 + y^2)) 

       val <- val1*p + (1 - p)*val2

       val[ u <= 0 | v <= 0 ] <- 0
       val[ u >= 1 | v >= 1] <- 0

      val
})

setMethod("dcopula","kimeldorf.sampson.copula",function(copula, u, v) {
   delta <- as.vector(copula@parameters[1])
   (1 + delta)*(u*v)^(-(delta + 1))*
       (u^(-delta) + v^(-delta) -1) ^(-2-1/delta)
})


setMethod("dcopula","husler.reiss.copula", function(copula, u, v) {
    d <- as.vector(copula@parameters[1])
    z <- log(u)/log(v)

    ff <- pnorm(1/d + 1/2*d*log(1/z)) *  pnorm(1/d + 1/2*d*log(z)) 

   1/u/v * pcopula(copula, u,v) * (ff + 1/2 * d/(-log(v))*dnorm(1/d + 1/2*d*log(z)))    
})


setMethod("dcopula","frank.copula", function(copula, u, v) {
     delta <- as.vector(copula@parameters[1])

     eta <- 1 - exp(-delta)
     delta *eta *exp( - delta *(u+v))/(eta 
        -(1 - exp(-delta*u))*(1 - exp(-delta*v)))^2   
})


setMethod("dcopula","bb4.copula", function(copula, u, v) {
       pu <- Phi(copula,u)
       pv <- Phi(copula,v)

       t <- pu/(pu+pv)
       ad0 <- A(copula, t)
       ad1 <- AfirstDer(copula, t)
       ad2 <- AsecondDer(copula, t)
       t2 <- ad0*(pu + pv)
       val <-  ((-(ad2*InvPhifirstDer(copula,t2)*pu*pv) + 
    InvPhisecondDer(copula,t2)*(pu+pv)*
       (ad0*(ad0 - ad1)*pu^2 + 
         (2*ad0^2 - ad1^2)*pu*pv + 
         ad0*(ad0 + ad1)*pv^2))*
    PhiDer(copula,u)*PhiDer(copula,v))/
  (pu+pv)^3

       val[is.na(val)] <- 0
       val
})



setMethod("dcopula","gumbel.copula",
     function(copula, u, v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    t <- log(x)/log(x * y)
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    ad2 <- AsecondDer(copula, t)
    lxy <- log(x * y)
    lx <- log(x)
    val <- ad0^2/(x * y) + (ad0 * ad1)/(x * y) + (ad2 * lx^2)/(x * y * lxy^3) - (ad2 * lx)/(x * y * lxy^2) +
        (ad1^2 * lx^2)/(x * y * lxy^2) - (2 * ad0 * ad1 * lx)/(x * y * lxy) - (ad1^2 * lx)/(x * y * lxy)
    val <- val * pcopula(copula, x, y)
    val[u == 0 | v == 0 | v == 1 | u == 1] <- 0
    #   val[is.na(val)] <- NA
    val
} )

setMethod("dcopula","galambos.copula",
     function(copula, u, v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    t <- log(x)/log(x * y)
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    ad2 <- AsecondDer(copula, t)
    lxy <- log(x * y)
    lx <- log(x)
    val <- ad0^2/(x * y) + (ad0 * ad1)/(x * y) + (ad2 * lx^2)/(x * y * lxy^3) - (ad2 * lx)/(x * y * lxy^2) +
        (ad1^2 * lx^2)/(x * y * lxy^2) - (2 * ad0 * ad1 * lx)/(x * y * lxy) - (ad1^2 * lx)/(x * y * lxy)
    val <- val * pcopula(copula, x, y)
    val[u == 0 | v == 0 | v == 1 | u == 1] <- 0
    #   val[is.na(val)] <- NA
    val
} )

setMethod("dcopula","tawn.copula",
     function(copula, u, v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    t <- log(x)/log(x * y)
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    ad2 <- AsecondDer(copula, t)
    lxy <- log(x * y)
    lx <- log(x)
    val <- ad0^2/(x * y) + (ad0 * ad1)/(x * y) + (ad2 * lx^2)/(x * y * lxy^3) - (ad2 * lx)/(x * y * lxy^2) +
        (ad1^2 * lx^2)/(x * y * lxy^2) - (2 * ad0 * ad1 * lx)/(x * y * lxy) - (ad1^2 * lx)/(x * y * lxy)
    val <- val * pcopula(copula, x, y)
    val[u == 0 | v == 0 | v == 1 | u == 1] <- 0
    #   val[is.na(val)] <- NA
    val
} )

setMethod("dcopula","bb5.copula",
     function(copula, u, v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    t <- log(x)/log(x * y)
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    ad2 <- AsecondDer(copula, t)
    lxy <- log(x * y)
    lx <- log(x)
    val <- ad0^2/(x * y) + (ad0 * ad1)/(x * y) + (ad2 * lx^2)/(x * y * lxy^3) - (ad2 * lx)/(x * y * lxy^2) +
        (ad1^2 * lx^2)/(x * y * lxy^2) - (2 * ad0 * ad1 * lx)/(x * y * lxy) - (ad1^2 * lx)/(x * y * lxy)
    val <- val * pcopula(copula, x, y)
    val[u == 0 | v == 0 | v == 1 | u == 1] <- 0
    #   val[is.na(val)] <- NA
    val
} )


setMethod("dcopula","bb1.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    val <- InvPhisecondDer(copula, val)
    val <- val * PhiDer(copula, u) * PhiDer(copula, v)
    val
} )

setMethod("dcopula","bb2.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    val <- InvPhisecondDer(copula, val)
    val <- val * PhiDer(copula, u) * PhiDer(copula, v)
    val
} )

setMethod("dcopula","bb3.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    val <- InvPhisecondDer(copula, val)
    val <- val * PhiDer(copula, u) * PhiDer(copula, v)
    val
} )

setMethod("dcopula","bb6.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    val <- InvPhisecondDer(copula, val)
    val <- val * PhiDer(copula, u) * PhiDer(copula, v)
    val
} )

setMethod("dcopula","bb7.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    val <- InvPhisecondDer(copula, val)
    val <- val * PhiDer(copula, u) * PhiDer(copula, v)
    val
} )

setMethod("dcopula","joe.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    val <- InvPhisecondDer(copula, val)
    val <- val * PhiDer(copula, u) * PhiDer(copula, v)
    val
} )


########################## 
### Functions pcopula  ### 


setMethod("pcopula","normal.copula", function(copula, u, v) {
       delta <- as.vector(copula@parameters[1])

      good <- u < 1 & u > 0 & v < 1 & v > 0      
      ug <- u[good]
      vg <- v[good]      
      if (sum(good) > 0) 
      {
      	tlen <- max(length(ug),length(vg))
       	x <- rep(qnorm(ug),length.out = tlen)
      	 y <- rep(qnorm(vg),length.out = tlen)           
       	XX <- matrix(ncol = 2, nrow= length(x), data = c(x,y), byrow = F)     
#       valg <- pmvnorm(XX, rho = delta)
       	ff <- function(x) p2Dnorm(x,rho=delta)
       	valg <- apply(XX,1,ff)    
       }
       val <- vector(length =   length(good), mode = "numeric")
       val[1:length(good)] <- NA    
       if (sum(good) > 0) { val[good] <- valg}
       val[ u <= 0 | v <= 0 ] <- 0
       val[ u >= 1 & v >= 1] <- 1
       val[is.na(val) & u >= 1] <- v[is.na(val) & u >= 1]
       val[is.na(val) & v >= 1] <- u[is.na(val) & v >= 1]
       val
} )


setMethod("pcopula","frank.copula", function(copula, u, v) {
   delta <- as.vector(copula@parameters[1])
   eta <- 1 - exp(-delta)
   -1/delta*log((eta - (1 - exp(-delta*u))*(1 - exp(-delta*v)))/eta)
})


setMethod("pcopula","kimeldorf.sampson.copula", function(copula, u, v) {
   delta <- as.vector(copula@parameters[1])
   (u^(-delta) + v^(-delta) -1) ^(-1/delta)
})

setMethod("pcopula","normal.mix.copula", function(copula, u, v) {
    p <- as.vector(copula@parameters[1])
    delta1 <- as.vector(copula@parameters[2])
   delta2 <- as.vector(copula@parameters[3])

      tlen <- max(length(u),length(v))
      u <- rep(u,length.out = tlen)
      v <- rep(v,length.out = tlen)

      good <- u < 1 & u > 0 & v < 1 & v > 0      
      ug <- u[good]        
      vg <- v[good]      

      x <- qnorm(ug)       
      y <- qnorm(vg)     

       XX <- matrix(ncol = 2, nrow= length(x), data = c(x,y), byrow = F)

     #  cat(ug,"\n")       
     #  cat(vg,"\n")       
     #  cat(sum(is.na(XX)),"\n")

#       valg1 <- pmvnorm(XX, rho = delta1)           
       ff <- function(x) p2Dnorm(x,rho=delta1)
       valg1 <- apply(XX,1,ff)    
#       valg2 <- pmvnorm(XX, rho = delta2)
       ff <- function(x) p2Dnorm(x,rho=delta2)
       valg2 <- apply(XX,1,ff)    

       valg <- p * valg1 + (1 -p) * valg2

       val <- vector(length = length(good), mode = "numeric")
       val[1:length(good)] <- NA    
       val[good] <- valg
       val[ u <= 0 | v <= 0 ] <- 0
       val[ u >= 1 & v >= 1] <- 1
       val[is.na(val) & u >= 1] <- v[is.na(val) & u >= 1]
       val[is.na(val) & v >= 1] <- u[is.na(val) & v >= 1]
       val
})


setMethod("pcopula","gumbel.copula",
     function(copula, u, v) {
    t <- log(u)/log(u * v)
    val <- exp(log(u * v) * A(copula, t))
    val[u <= 0 | v <= 0] <- 0
    val[u >= 1 & v >= 1] <- 1
    val[is.na(val) & u >= 1] <- v[is.na(val) & u >= 1]
    val[is.na(val) & v >= 1] <- u[is.na(val) & v >= 1]
    val
} )

setMethod("pcopula","galambos.copula",
     function(copula, u, v) {
    t <- log(u)/log(u * v)
    val <- exp(log(u * v) * A(copula, t))
    val[u <= 0 | v <= 0] <- 0
    val[u >= 1 & v >= 1] <- 1
    val[is.na(val) & u >= 1] <- v[is.na(val) & u >= 1]
    val[is.na(val) & v >= 1] <- u[is.na(val) & v >= 1]
    val
} )


setMethod("pcopula","husler.reiss.copula",function(copula, u, v) {
  d <- as.vector(copula@parameters[1])
  uu <- -log(u)
  vv <- -log(v)
  ft1 <- pnorm(1/d + 1/2*d*log(uu/vv))
  ft2 <- pnorm(1/d + 1/2*d*log(vv/uu))

  exp(-uu*ft1 - vv *ft2)
 })

setMethod("pcopula","tawn.copula",
     function(copula, u, v) {
    t <- log(u)/log(u * v)
    val <- exp(log(u * v) * A(copula, t))
    val[u <= 0 | v <= 0] <- 0
    val[u >= 1 & v >= 1] <- 1
    val[is.na(val) & u >= 1] <- v[is.na(val) & u >= 1]
    val[is.na(val) & v >= 1] <- u[is.na(val) & v >= 1]
    val
} )

setMethod("pcopula","bb5.copula",
     function(copula, u, v) {
    t <- log(u)/log(u * v)
    val <- exp(log(u * v) * A(copula, t))
    val[u <= 0 | v <= 0] <- 0
    val[u >= 1 & v >= 1] <- 1
    val[is.na(val) & u >= 1] <- v[is.na(val) & u >= 1]
    val[is.na(val) & v >= 1] <- u[is.na(val) & v >= 1]
    val
} )

setMethod("pcopula","bb1.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    InvPhi(copula, val)
} )

setMethod("pcopula","bb2.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    InvPhi(copula, val)
} )

setMethod("pcopula","bb3.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    InvPhi(copula, val)
} )

setMethod("pcopula","bb6.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    InvPhi(copula, val)
} )

setMethod("pcopula","bb7.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    InvPhi(copula, val)
} )

setMethod("pcopula","joe.copula",
     function(copula, u, v) {
    val <- Phi(copula, u) + Phi(copula, v)
    InvPhi(copula, val)
} )

setMethod("pcopula","bb4.copula",function(copula, u,v) {
      pu <- Phi(copula,u)
      pv <- Phi(copula,v)
      val <- A(copula,pu/(pu+pv))
     InvPhi(copula,(pu+pv)*val)
})


setMethod("pcopula","empirical.copula", function(copula, u, v) 
{
    data <- copula
    obs.x <- c(0,data@x,1)
    obs.y <- c(0,data@y,1)
    slistX <- sort.list(obs.x)
   slistY <- sort.list(obs.y)
   xsort <- as.double(obs.x[slistX])
   ysort <- as.double(obs.y[slistY])
   ind <- as.integer(match(slistX,slistY) - 1)
   nsamp <- as.integer(length(xsort))

   N.in <- length(u)
    good.in.values <- (u >= 0) & (u <= 1) & (v >= 0) & (v <= 1)
    val <- vector (length = N.in, mode = "numeric")
    val[!good.in.values] <- NA
    xx <- as.double(u[good.in.values])
    yy <- as.double(v[good.in.values]) 
    N <- as.integer(length(xx))
    val1 <- as.double(rep(-1,N))

   outC <- .C("interpcopR", xsort,ysort,ind, nsamp,xx,yy,
                N,val1)
   val[good.in.values] <- outC[[8]]
    val
})


########################## 
### Functions rcopula  ### 


setMethod("rcopula", "normal.copula" , function(copula, n) {
       delta <-  copula@parameters[1]
       bvsn <- rmvnorm(n,cov=matrix(c(1,delta,delta,1),byrow=T,ncol=2))
       x <- pnorm(bvsn[,1])
       y <- pnorm(bvsn[,2])
       sim <- list(x = x, y = y)
       sim
})

setMethod("rcopula","frank.copula", function(copula, n) {
   # inverse EXISTS!!
   Usamp <- runif(n)
   Vsamp <- CondCinverse.b3.copula(copula,runif(n),Usamp)
   val <- list(x = Usamp, y = Vsamp)
   val
})

setMethod("rcopula","kimeldorf.sampson.copula", function(copula, n) {
   # inverse EXISTS!!
   Usamp <- runif(n)
   Vsamp <- CondCinverse(copula,runif(n),Usamp)
   val <- list(x = Usamp, y = Vsamp)
   val
})

setMethod("rcopula","gumbel.copula",
     function(copula, n) {
    # generate z from distribution h
    # using rejection method
    zrand <- vector(length = n, mode = "numeric")
    xgenerated <- vector(length = n, mode = "logical")
    lg <- n
    #set c to 1/2 true for symmetric copula
    cc <- Hderiv(copula, 1/2)
    stopnow <- F
    while(!stopnow) {
        usamp <- runif(lg)
        ysamp <- runif(lg)
        ykeep <- (usamp <= Hderiv(copula, ysamp)/cc)
        #indexes of xrand vector to store "good" generated values
        toind <- (c(1:n)[!xgenerated])[ykeep]
        zrand[toind] <- ysamp[ykeep]
        xgenerated[toind] <- T
        #number of rvs left to generate:
        lg <- n - sum(xgenerated)
        if(lg == 0)
            stopnow <- T
    }
    z <- zrand
    u <- runif(n)
    # cat(length(z), length(AsecondDer(copula,z)))
    pz <- (z * (1 - z) * AsecondDer(copula, z))/Hderiv(copula, z)/A(copula, z)
    w <- NULL
    w[1:n] <- 0
    nn <- sum(u <= pz)
    w[u <= pz] <- runif(nn)
    w[u > pz] <- runif(n - nn) * runif(n - nn)
    xx <- exp((z * log(w))/A(copula, z))
    yy <- exp(((1 - z) * log(w))/A(copula, z))
    val <- list(x = xx, y = yy)
    val
} )

setMethod("rcopula","galambos.copula",
     function(copula, n) {
    # generate z from distribution h
    # using rejection method
    zrand <- vector(length = n, mode = "numeric")
    xgenerated <- vector(length = n, mode = "logical")
    lg <- n
    #set c to 1/2 true for symmetric copula
    cc <- Hderiv(copula, 1/2)
    stopnow <- F
    while(!stopnow) {
        usamp <- runif(lg)
        ysamp <- runif(lg)
        ykeep <- (usamp <= Hderiv(copula, ysamp)/cc)
        #indexes of xrand vector to store "good" generated values
        toind <- (c(1:n)[!xgenerated])[ykeep]
        zrand[toind] <- ysamp[ykeep]
        xgenerated[toind] <- T
        #number of rvs left to generate:
        lg <- n - sum(xgenerated)
        if(lg == 0)
            stopnow <- T
    }
    z <- zrand
    u <- runif(n)
    # cat(length(z), length(AsecondDer(copula,z)))
    pz <- (z * (1 - z) * AsecondDer(copula, z))/Hderiv(copula, z)/A(copula, z)
    w <- NULL
    w[1:n] <- 0
    nn <- sum(u <= pz)
    w[u <= pz] <- runif(nn)
    w[u > pz] <- runif(n - nn) * runif(n - nn)
    xx <- exp((z * log(w))/A(copula, z))
    yy <- exp(((1 - z) * log(w))/A(copula, z))
    val <- list(x = xx, y = yy)
    val
} )

setMethod("rcopula","husler.reiss.copula",
     function(copula, n) {
    # generate z from distribution h
    # using rejection method
    zrand <- vector(length = n, mode = "numeric")
    xgenerated <- vector(length = n, mode = "logical")
    lg <- n
    #set c to 1/2 true for symmetric copula
    cc <- Hderiv(copula, 1/2)
    stopnow <- F
    while(!stopnow) {
        usamp <- runif(lg)
        ysamp <- runif(lg)
        ykeep <- (usamp <= Hderiv(copula, ysamp)/cc)
        #indexes of xrand vector to store "good" generated values
        toind <- (c(1:n)[!xgenerated])[ykeep]
        zrand[toind] <- ysamp[ykeep]
        xgenerated[toind] <- T
        #number of rvs left to generate:
        lg <- n - sum(xgenerated)
        if(lg == 0)
            stopnow <- T
    }
    z <- zrand
    u <- runif(n)
    # cat(length(z), length(AsecondDer(copula,z)))
    pz <- (z * (1 - z) * AsecondDer(copula, z))/Hderiv(copula, z)/A(copula, z)
    w <- NULL
    w[1:n] <- 0
    nn <- sum(u <= pz)
    w[u <= pz] <- runif(nn)
    w[u > pz] <- runif(n - nn) * runif(n - nn)
    xx <- exp((z * log(w))/A(copula, z))
    yy <- exp(((1 - z) * log(w))/A(copula, z))
    val <- list(x = xx, y = yy)
    val
} )

setMethod("rcopula","tawn.copula",
     function(copula, n) {
    # generate z from distribution h
    # using rejection method
    zrand <- vector(length = n, mode = "numeric")
    xgenerated <- vector(length = n, mode = "logical")
    lg <- n
    #set c to 1/2 true for symmetric copula
    cc <- Hderiv(copula, 1/2)
    stopnow <- F
    while(!stopnow) {
        usamp <- runif(lg)
        ysamp <- runif(lg)
        ykeep <- (usamp <= Hderiv(copula, ysamp)/cc)
        #indexes of xrand vector to store "good" generated values
        toind <- (c(1:n)[!xgenerated])[ykeep]
        zrand[toind] <- ysamp[ykeep]
        xgenerated[toind] <- T
        #number of rvs left to generate:
        lg <- n - sum(xgenerated)
        if(lg == 0)
            stopnow <- T
    }
    z <- zrand
    u <- runif(n)
    # cat(length(z), length(AsecondDer(copula,z)))
    pz <- (z * (1 - z) * AsecondDer(copula, z))/Hderiv(copula, z)/A(copula, z)
    w <- NULL
    w[1:n] <- 0
    nn <- sum(u <= pz)
    w[u <= pz] <- runif(nn)
    w[u > pz] <- runif(n - nn) * runif(n - nn)
    xx <- exp((z * log(w))/A(copula, z))
    yy <- exp(((1 - z) * log(w))/A(copula, z))
    val <- list(x = xx, y = yy)
    val
} )

setMethod("rcopula","bb5.copula",
     function(copula, n) {
    # generate z from distribution h
    # using rejection method
    zrand <- vector(length = n, mode = "numeric")
    xgenerated <- vector(length = n, mode = "logical")
    lg <- n
    #set c to 1/2 true for symmetric copula
    cc <- Hderiv(copula, 1/2)
    stopnow <- F
    while(!stopnow) {
        usamp <- runif(lg)
        ysamp <- runif(lg)
        ykeep <- (usamp <= Hderiv(copula, ysamp)/cc)
        #indexes of xrand vector to store "good" generated values
        toind <- (c(1:n)[!xgenerated])[ykeep]
        zrand[toind] <- ysamp[ykeep]
        xgenerated[toind] <- T
        #number of rvs left to generate:
        lg <- n - sum(xgenerated)
        if(lg == 0)
            stopnow <- T
    }
    z <- zrand
    u <- runif(n)
    # cat(length(z), length(AsecondDer(copula,z)))
    pz <- (z * (1 - z) * AsecondDer(copula, z))/Hderiv(copula, z)/A(copula, z)
    w <- NULL
    w[1:n] <- 0
    nn <- sum(u <= pz)
    w[u <= pz] <- runif(nn)
    w[u > pz] <- runif(n - nn) * runif(n - nn)
    xx <- exp((z * log(w))/A(copula, z))
    yy <- exp(((1 - z) * log(w))/A(copula, z))
    val <- list(x = xx, y = yy)
    val
} )


setMethod("rcopula","bb1.copula",
     function(copula, n) {
    #figure out family:
    famstr <- substring(class(copula)[1], 1, 3)
    implem.fam <- c("bb1", "bb2", "bb3", "bb6", "bb7")
    famID <- c(1, 2, 3, 6, 7)
    ind <- charmatch(famstr, implem.fam)
    family <- as.integer(famID[ind])
    #use Genest MacKay 1996algorithm
    u <- runif(n)
    t <- runif(n)
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
    # use C code to compute the inverse of secon deriv. of Phi
    x <- PhiDer(copula, u)/t
    # this is temporary!!!
    infx <- is.inf(x)
    while(sum(infx) > 0) {
        warning(paste(" The derivative of Phi is Infinity in ", sum(infx), " points",
            " Simulations results maybe unreliable", sep = ""))
        tinf <- sum(infx)
        x[infx] <- PhiDer(copula, runif(tinf)/runif(tinf))
        infx <- is.inf(x)
    }
    Xtemp <- as.double(sort(x))
    result <- vector(length = n)
    result[1:n] <- as.double(-999)
    epsX <- as.double(1e-010)
    epsY <- as.double(1e-010)
    conv <- as.integer(rep(-9, n))
    maxit <- as.integer(100)
    cout <- .C("inverseDerivPhiBB",
        Xtemp,
        as.integer(n),
        delta,
        theta,
        result,
        family,
        epsX,
        epsY,
        maxit,
        conv)
    result <- cout[[5]]
    conv <- cout[[10]]
    if(sum(conv <= 0) != 0)
        warning("Inverting derivative of Phi is unsuccessful. Simulation results might be unreliable")
    w <- result
    w[sort.list(x)] <- result
    Vtemp <- Phi(copula, w) - Phi(copula, u)
    v <- InvPhi(copula, Vtemp)
    val <- list(x = u, y = v)
    val
} )

setMethod("rcopula","bb2.copula",
     function(copula, n) {
    #figure out family:
    famstr <- substring(class(copula)[1], 1, 3)
    implem.fam <- c("bb1", "bb2", "bb3", "bb6", "bb7")
    famID <- c(1, 2, 3, 6, 7)
    ind <- charmatch(famstr, implem.fam)
    family <- as.integer(famID[ind])
    #use Genest MacKay 1996algorithm
    u <- runif(n)
    t <- runif(n)
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
    # use C code to compute the inverse of secon deriv. of Phi
    x <- PhiDer(copula, u)/t
    # this is temporary!!!
    infx <- is.inf(x)
    while(sum(infx) > 0) {
        warning(paste(" The derivative of Phi is Infinity in ", sum(infx), " points",
            " Simulations results maybe unreliable", sep = ""))
        tinf <- sum(infx)
        x[infx] <- PhiDer(copula, runif(tinf)/runif(tinf))
        infx <- is.inf(x)
    }
    Xtemp <- as.double(sort(x))
    result <- vector(length = n)
    result[1:n] <- as.double(-999)
    epsX <- as.double(1e-010)
    epsY <- as.double(1e-010)
    conv <- as.integer(rep(-9, n))
    maxit <- as.integer(100)
    cout <- .C("inverseDerivPhiBB",
        Xtemp,
        as.integer(n),
        delta,
        theta,
        result,
        family,
        epsX,
        epsY,
        maxit,
        conv)
    result <- cout[[5]]
    conv <- cout[[10]]
    if(sum(conv <= 0) != 0)
        warning("Inverting derivative of Phi is unsuccessful. Simulation results might be unreliable")
    w <- result
    w[sort.list(x)] <- result
    Vtemp <- Phi(copula, w) - Phi(copula, u)
    v <- InvPhi(copula, Vtemp)
    val <- list(x = u, y = v)
    val
} )

setMethod("rcopula","bb3.copula",
     function(copula, n) {
    #figure out family:
    famstr <- substring(class(copula)[1], 1, 3)
    implem.fam <- c("bb1", "bb2", "bb3", "bb6", "bb7")
    famID <- c(1, 2, 3, 6, 7)
    ind <- charmatch(famstr, implem.fam)
    family <- as.integer(famID[ind])
    #use Genest MacKay 1996algorithm
    u <- runif(n)
    t <- runif(n)
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
    # use C code to compute the inverse of secon deriv. of Phi
    x <- PhiDer(copula, u)/t
    # this is temporary!!!
    infx <- is.inf(x)
    while(sum(infx) > 0) {
        warning(paste(" The derivative of Phi is Infinity in ", sum(infx), " points",
            " Simulations results maybe unreliable", sep = ""))
        tinf <- sum(infx)
        x[infx] <- PhiDer(copula, runif(tinf)/runif(tinf))
        infx <- is.inf(x)
    }
    Xtemp <- as.double(sort(x))
    result <- vector(length = n)
    result[1:n] <- as.double(-999)
    epsX <- as.double(1e-010)
    epsY <- as.double(1e-010)
    conv <- as.integer(rep(-9, n))
    maxit <- as.integer(100)
    cout <- .C("inverseDerivPhiBB",
        Xtemp,
        as.integer(n),
        delta,
        theta,
        result,
        family,
        epsX,
        epsY,
        maxit,
        conv)
    result <- cout[[5]]
    conv <- cout[[10]]
    if(sum(conv <= 0) != 0)
        warning("Inverting derivative of Phi is unsuccessful. Simulation results might be unreliable")
    w <- result
    w[sort.list(x)] <- result
    Vtemp <- Phi(copula, w) - Phi(copula, u)
    v <- InvPhi(copula, Vtemp)
    val <- list(x = u, y = v)
    val
} )

setMethod("rcopula","bb6.copula",
     function(copula, n) {
    #figure out family:
    famstr <- substring(class(copula)[1], 1, 3)
    implem.fam <- c("bb1", "bb2", "bb3", "bb6", "bb7")
    famID <- c(1, 2, 3, 6, 7)
    ind <- charmatch(famstr, implem.fam)
    family <- as.integer(famID[ind])
    #use Genest MacKay 1996algorithm
    u <- runif(n)
    t <- runif(n)
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
    # use C code to compute the inverse of secon deriv. of Phi
    x <- PhiDer(copula, u)/t
    # this is temporary!!!
    infx <- is.inf(x)
    while(sum(infx) > 0) {
        warning(paste(" The derivative of Phi is Infinity in ", sum(infx), " points",
            " Simulations results maybe unreliable", sep = ""))
        tinf <- sum(infx)
        x[infx] <- PhiDer(copula, runif(tinf)/runif(tinf))
        infx <- is.inf(x)
    }
    Xtemp <- as.double(sort(x))
    result <- vector(length = n)
    result[1:n] <- as.double(-999)
    epsX <- as.double(1e-010)
    epsY <- as.double(1e-010)
    conv <- as.integer(rep(-9, n))
    maxit <- as.integer(100)
    cout <- .C("inverseDerivPhiBB",
        Xtemp,
        as.integer(n),
        delta,
        theta,
        result,
        family,
        epsX,
        epsY,
        maxit,
        conv)
    result <- cout[[5]]
    conv <- cout[[10]]
    if(sum(conv <= 0) != 0)
        warning("Inverting derivative of Phi is unsuccessful. Simulation results might be unreliable")
    w <- result
    w[sort.list(x)] <- result
    Vtemp <- Phi(copula, w) - Phi(copula, u)
    v <- InvPhi(copula, Vtemp)
    val <- list(x = u, y = v)
    val
} )

setMethod("rcopula","bb7.copula",
     function(copula, n) {
    #figure out family:
    famstr <- substring(class(copula)[1], 1, 3)
    implem.fam <- c("bb1", "bb2", "bb3", "bb6", "bb7")
    famID <- c(1, 2, 3, 6, 7)
    ind <- charmatch(famstr, implem.fam)
    family <- as.integer(famID[ind])
    #use Genest MacKay 1996algorithm
    u <- runif(n)
    t <- runif(n)
    theta <- as.vector(copula@parameters[1])
    delta <- as.vector(copula@parameters[2])
    # use C code to compute the inverse of secon deriv. of Phi
    x <- PhiDer(copula, u)/t
    # this is temporary!!!
    infx <- is.inf(x)
    while(sum(infx) > 0) {
        warning(paste(" The derivative of Phi is Infinity in ", sum(infx), " points",
            " Simulations results maybe unreliable", sep = ""))
        tinf <- sum(infx)
        x[infx] <- PhiDer(copula, runif(tinf)/runif(tinf))
        infx <- is.inf(x)
    }
    Xtemp <- as.double(sort(x))
    result <- vector(length = n)
    result[1:n] <- as.double(-999)
    epsX <- as.double(1e-010)
    epsY <- as.double(1e-010)
    conv <- as.integer(rep(-9, n))
    maxit <- as.integer(100)
    cout <- .C("inverseDerivPhiBB",
        Xtemp,
        as.integer(n),
        delta,
        theta,
        result,
        family,
        epsX,
        epsY,
        maxit,
        conv)
    result <- cout[[5]]
    conv <- cout[[10]]
    if(sum(conv <= 0) != 0)
        warning("Inverting derivative of Phi is unsuccessful. Simulation results might be unreliable")
    w <- result
    w[sort.list(x)] <- result
    Vtemp <- Phi(copula, w) - Phi(copula, u)
    v <- InvPhi(copula, Vtemp)
    val <- list(x = u, y = v)
    val
} )

setMethod("rcopula","joe.copula", function(copula, n) {
       theta <- as.vector(copula@parameters[1])
       rcopula(bb6.copula(theta,1.0), n)
    })


setMethod("rcopula","bb4.copula", function(copula, n) {
   # generate z from distribution h
   # using rejection method

   zrand <- vector(length = n, mode = "numeric")
   xgenerated <- vector(length = n, mode = "logical")
   lg <- n
   #set c to 1/2 true for symmetric copula
   cc <- Hderiv(copula,1/2)
   stopnow <- F
   while(!stopnow) {
     usamp <- runif(lg)
     ysamp <- runif(lg)
     ykeep <- (usamp <= Hderiv(copula,ysamp)/cc)
     #indexes of xrand vector to store "good" generated values
     toind <- (c(1:n)[!xgenerated])[ykeep]
     zrand[toind] <- ysamp[ykeep]
     xgenerated[toind] <- T
     #number of rvs left to generate:
     lg <- n - sum(xgenerated)
     if (lg == 0)  stopnow <- T
   }

   z <- zrand 
   u <- runif(n)

 # cat(length(z), length(AsecondDer(copula,z)))
   pz <- z*(1-z)*AsecondDer(copula,z)/Hderiv(copula,z)/A(copula,z)

   w <- NULL
   w[1:n] <- 0
   nn <- sum(u <=  pz)
   w[u <=  pz] <- runif(nn)

   # we need to generate n - nn r.v from 
   # archimedian copula with phi
      PHIU <- runif(n - nn)
      PHIT <- runif(n - nn)
      PHIW <- DerPhi.minus1(copula, PhiDer(copula,PHIU)/PHIT)

   w[u > pz] <- PHIW

   pw <- Phi(copula,w)
   xx <- InvPhi(copula, z*pw/A(copula,z))
   yy <- InvPhi(copula, (1 - z)*pw/A(copula,z))

   val <- list(x = xx, v = yy)
   val   
})


setMethod("rcopula","normal.mix.copula", function(copula, n) {
          p <- as.vector(copula@parameters[1])
    delta1 <- as.vector(copula@parameters[2])
   delta2 <- as.vector(copula@parameters[3])

       bvsn1 <- rmvnorm(n, d = 2, rho = delta1)
       bvsn2 <- rmvnorm(n, d = 2, rho = delta2)

       uu <- runif(n)

      bvsn <- bvsn1
      bvsn[uu > p,] <- bvsn2[uu > p,]      

       x <- pnorm(bvsn[,1])
       y <- pnorm(bvsn[,2])
       sim <- list(x = x, y = y)
       sim  
})



######################### 
### Functions Hderiv  ### 


setMethod("Hderiv","galambos.copula",
     function(copula, z) {
    ad0 <- A(copula, z)
    ad1 <- AfirstDer(copula, z)
    ad2 <- AsecondDer(copula, z)
    1 + ((1 - 2 * z) * ad1)/ad0 + z * (1 - z) * (ad2/ad0 - ad1^2/ad0^2)
} )

setMethod("Hderiv","husler.reiss.copula",
     function(copula, z) {
    ad0 <- A(copula, z)
    ad1 <- AfirstDer(copula, z)
    ad2 <- AsecondDer(copula, z)
    1 + ((1 - 2 * z) * ad1)/ad0 + z * (1 - z) * (ad2/ad0 - ad1^2/ad0^2)
} )

setMethod("Hderiv","bb5.copula",
     function(copula, z) {
    ad0 <- A(copula, z)
    ad1 <- AfirstDer(copula, z)
    ad2 <- AsecondDer(copula, z)
    1 + ((1 - 2 * z) * ad1)/ad0 + z * (1 - z) * (ad2/ad0 - ad1^2/ad0^2)
} )




############################### 
### Functions Kendalls.tau  ### 


setMethod("Kendalls.tau","gumbel.copula", function(copula, tol = 1e-5)
 {
    f <- function(t, copula)
    {
        (t * (1 - t) * AsecondDer(copula, t))/A(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    t1
} )

setMethod("Kendalls.tau","galambos.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        (t * (1 - t) * AsecondDer(copula, t))/A(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    t1
} )

setMethod("Kendalls.tau","husler.reiss.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        (t * (1 - t) * AsecondDer(copula, t))/A(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    t1
} )

setMethod("Kendalls.tau","tawn.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        (t * (1 - t) * AsecondDer(copula, t))/A(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    t1
} )

setMethod("Kendalls.tau","bb5.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        (t * (1 - t) * AsecondDer(copula, t))/A(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    t1
} )

setMethod("Kendalls.tau","normal.copula", function(copula, tol = 1e-5) {
	  delta <- as.vector(copula@parameters[1])
     return(2*asin(delta)/pi)	
})

setMethod("Kendalls.tau","kimeldorf.sampson.copula", function(copula, tol = 1e-5) {
	  delta <- as.vector(copula@parameters[1])
     return(delta/(delta+2))	
})

setMethod("Kendalls.tau","bb1.copula", function(copula, tol = 1e-5) {
	  theta <- as.vector(copula@parameters[1])
	  delta <- as.vector(copula@parameters[2])
	  return(1 - 2/delta/(theta + 2))	
})

setMethod("Kendalls.tau","bb6.copula", function(copula, tol = 1e-5) {
	  delta <- as.vector(copula@parameters[1])
     return((delta-1)/(delta))	
})

setMethod("Kendalls.tau","bb4.copula", function(copula, tol = 1e-5){
	  theta <- as.vector(copula@parameters[1])
	  delta <- as.vector(copula@parameters[2])
	  tauphi <- 1 - 4/(4 + 2*theta)
	   
	   f <- function(t, copula) 
	   {	   
		    t*(1 - t)*AsecondDer(copula,t)/A(copula,t)	
	}   		
	   t1 <- integrate(f,subdivisions=100, lower = 0.0, 
                        upper = 1.0, copula = copula)$value      
	   return(t1+(1-t1)*tauphi)	
})

setMethod("Kendalls.tau","frank.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        Phi(copula, t)/PhiDer(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    1 + 4 * t1
} )

setMethod("Kendalls.tau","bb2.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        Phi(copula, t)/PhiDer(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    1 + 4 * t1
} )

setMethod("Kendalls.tau","bb3.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        Phi(copula, t)/PhiDer(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    1 + 4 * t1
} )

setMethod("Kendalls.tau","bb6.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        Phi(copula, t)/PhiDer(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    1 + 4 * t1
} )

setMethod("Kendalls.tau","bb7.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        Phi(copula, t)/PhiDer(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    1 + 4 * t1
} )

setMethod("Kendalls.tau","joe.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        Phi(copula, t)/PhiDer(copula, t)
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    1 + 4 * t1
} )

setMethod("Kendalls.tau","normal.mix.copula", function(copula, tol = 1e-2) 
{          
		   assign("tempCop", copula, frame = 1)		
		   Cv <- function(x,y) {
			        x1 <- pcopula(tempCop,x,y)
              	 x2 <- dcopula(tempCop,x,y)
                  return(x1*x2)  
	      }
	  fun <- function(t, integrand, tol)
           { 
           	 val <- unlist(lapply(t, function(t, integrand, tol)
                 { integrate(f = integrand,
                    lower = 0, upper = 1, rel.tol = tol,
                    y = t)$value }
                , integrand, tol = tol))
               val
            }      
		   t1 <- integrate(f = fun,  subdivisions=100, lower = 0, rel.tol = tol,
                        upper = 1, integrand = Cv, tol = tol )$value
          val2 <- 4*t1 - 1
	    val2	
})


################################ 
### Functions Spearmans.rho  ### 


setMethod("Spearmans.rho","gumbel.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        1/(A(copula, t) + 1)^2
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","galambos.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        1/(A(copula, t) + 1)^2
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","husler.reiss.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        1/(A(copula, t) + 1)^2
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","tawn.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        1/(A(copula, t) + 1)^2
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","bb5.copula",
     function(copula, tol = 1e-5) {
    f <- function(t, copula)
    {
        1/(A(copula, t) + 1)^2
    }
    t1 <- integrate(f, subdivisions = 100, lower = 0., upper = 1., copula = copula)$value
    12 * t1 - 3
} )


setMethod("Spearmans.rho","normal.copula", function(copula, tol = 1e-5) {
	  delta <- as.vector(copula$parameters[1])
     return(6*asin(delta/2.0)/pi)	
})

setMethod("Spearmans.rho","normal.mix.copula", function(copula, tol = 1e-5) 
{    
		   assign("tempCop", copula, frame = 1)
		
	 Cv <- function(x,y) {
			        x1 <- pcopula(tempCop,x,y)
			        return(x1)  
	      }	 
	 fun <- function(t, integrand)
           {  val <- unlist(lapply(t, function(t, integrand)
                 { integrate(f = integrand,
                    lower = 0, upper = 1, 
                    y = t)$value }
                , integrand))
               val
            } 

	 t1 <- integrate(f = fun,  subdivisions=100, lower = 0, upper = 1, integrand = Cv )$value
        12*t1 - 3	
})

setMethod("Spearmans.rho","bb4.copula", function(copula, tol = 1e-5) 
{    
		   assign("tempCop", copula, frame = 1)
		
	 Cv <- function(x,y) {
			        x1 <- pcopula(tempCop,x,y)
			        return(x1)  
	      }

		 
	 fun <- function(t, integrand)
           {  val <- unlist(lapply(t, function(t, integrand)
                 { integrate(f = integrand,
                    lower = 0, upper = 1, 
                    y = t)$value }
                , integrand))
               val
            } 
	 t1 <- integrate(f = fun,  subdivisions=100, lower = 0, upper = 1, integrand = Cv )$value
          12*t1 - 3	
})


setMethod("Spearmans.rho","frank.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","kimeldorf.sampson.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","bb1.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","bb2.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","bb3.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","bb6.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","bb7.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )

setMethod("Spearmans.rho","joe.copula",
     function(copula, tol = 1e-5) {
    assign("tempCop", copula, frame = 1)
    Cv <- function(x, y)
    {
        x1 <- pcopula(tempCop, x, y)
        return(x1)
    }
    fun <- function(t, integrand)
    {
        val <- unlist(lapply(t, function(t, integrand)
        {
            integrate(f = integrand, lower = 0, upper = 1, y = t)$value
        }
        , integrand))
        val
    }
    t1 <- integrate(f = fun, subdivisions = 100, lower = 0, upper = 1, integrand = Cv)$value
    12 * t1 - 3
} )



############################# 
### Functions tail.index  ### 


setMethod("tail.index","gumbel.copula",
     function(copula, ...) {
    LI <- 0
    UI <- 2 - 2 * A(copula, 0.5)
    val <- c(LI, UI)
    names(val) <- c("lower.tail", "upper.tail")
    val
} )

setMethod("tail.index","galambos.copula",
     function(copula, ...) {
    LI <- 0
    UI <- 2 - 2 * A(copula, 0.5)
    val <- c(LI, UI)
    names(val) <- c("lower.tail", "upper.tail")
    val
} )

setMethod("tail.index","husler.reiss.copula",
     function(copula, ...) {
    LI <- 0
    UI <- 2 - 2 * A(copula, 0.5)
    val <- c(LI, UI)
    names(val) <- c("lower.tail", "upper.tail")
    val
} )

setMethod("tail.index","tawn.copula",
     function(copula, ...) {
    LI <- 0
    UI <- 2 - 2 * A(copula, 0.5)
    val <- c(LI, UI)
    names(val) <- c("lower.tail", "upper.tail")
    val
} )

setMethod("tail.index","bb5.copula",
     function(copula, ...) {
    LI <- 0
    UI <- 2 - 2 * A(copula, 0.5)
    val <- c(LI, UI)
    names(val) <- c("lower.tail", "upper.tail")
    val
} )




############################################################
# Plotting Functions for Copulas and Bivariate Distributions

persp.copula <- function(cop, n = 20, ...) {

     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
    divis <- divis[2:(n+1)]
        xmat <- rep(divis, n )
    ymat <- rep(divis, each = n)
    zmat <- dcopula(cop,  xmat,  ymat)

    val <- list(x = divis, y = divis, z = matrix(ncol = n, nrow = n, byrow = F, 
       data = zmat)   )

persp3d(divis, divis, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "u", ylab = "v", zlab = "Density")

#    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat),
#             z = as.vector(zmat)))    
#    guiPlot("Data Grid Surface", DataSetValues = val2, ...) 

    invisible(val)
}


persp.dcopula <- function(cop, n = 20, ...) {

     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
    divis <- divis[2:(n+1)]
        xmat <- rep(divis, n )
    ymat <- rep(divis, each = n)
    if (!inherits(cop,"empirical.copula")) { zmat <- dcopula(cop,  xmat,  ymat) }
    else {  bigobj <- kde2d(cop@x,cop@y,n = n, lims = c(range(divis),range(divis)))
            zmat <-   bigobj$z
    }

    val <- list(x = divis, y = divis, z = matrix(ncol = n, nrow = n, byrow = F, 
       data = zmat)   )

persp3d(divis, divis, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "Density")

#    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat),
#             z = as.vector(zmat)))    
#    guiPlot("Data Grid Surface", DataSetValues = val2, ...) 
    invisible(val)
}

persp.pcopula <- function(cop, n = 20, ...) {

     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
    divis <- divis[2:(n+1)]
        xmat <- rep(divis, n )
    ymat <- rep(divis, each = n)
    zmat <- pcopula(cop,  xmat,  ymat)

    val <- list(x = divis, y = divis, z = matrix(ncol = n, nrow = n, byrow = F, 
       data = zmat)   )

persp3d(divis, divis, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "Copula")

#    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat),
#             z = as.vector(zmat)))
#    guiPlot("Data Grid Surface", DataSetValues = val2, ...) 
    invisible(val)
}



contour.dcopula <- function(x,  ...) 
{args <- list(...)
  if(!"n"%in%names(args)) args$n <- 50
      n <- args$n
     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
    divis <- divis[2:(n+1)]
        xmat <- rep(divis, n )
    ymat <- rep(divis, each = n)
    zmat <- dcopula(x,  xmat,  ymat)

    val <- list(x = divis, y = divis, z = matrix(ncol = n, nrow = n, byrow = F,  data = zmat)   )
    contour(val)
# contour(x,y,val$z)
title("Contour Plot of the Copula Density", xlab="u",ylab="v")
    invisible(val)
}

contour.pcopula <- function(x, ...)
 { args <- list(...)
  if(!"n"%in%names(args)) args$n <- 50
      n <- args$n
     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
    divis <- divis[2:(n+1)]
        xmat <- rep(divis, n )
    ymat <- rep(divis, each = n)
    zmat <- pcopula(x,  xmat,  ymat)

    val <- list(x = divis, y = divis, z = matrix(ncol = n, nrow = n, byrow = F, data = zmat)   )
contour(val)
#	contour(x,y,val$z)
title("Contour Plot of the Copula", xlab="u",ylab="v")
    invisible(val)
}


############################### 
### Functions contour.plot  ###

# if(isGeneric("contour.plot"))removeGeneric("contour.plot")
# setGeneric("contour.plot",function(object, n = 100, nlevels = 10, add = F,  ...) standardGeneric("contour.plot"))


setMethod("contour.plot","normal.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","frank.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","kimeldorf.sampson.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","gumbel.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","galambos.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","husler.reiss.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","tawn.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb5.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","normal.mix.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb1.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb2.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb3.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb6.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb7.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","joe.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    # par(pin = c(6,6))
    contour(val, nlevels = nlevels, add = add, ..., xlab = "U", ylab = "V")
    invisible(val)
} )

setMethod("contour.plot","bb4.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
    contour(val, nlevels = nlevels, add = add, ..., xlab = "u", ylab = "v")
    invisible(val)
} )

setMethod("contour.plot","empirical.copula",
     function(object, n = 100, nlevels = 10, add = F,  ...){
    # oldpar <- par()
    # on.exit(par(oldpar))
    # cat(class(x))
    val <- list(x = (c(0:n)/n), y = (c(0:n)/n), z = matrix(ncol = n + 1, nrow = n + 1, data = NA))
    xmat <- rep(c(0:n)/n, n + 1)
    ymat <- rep(c(0:n)/n, each = n + 1)
    # ymat <- matrix(ncol = n+1, nrow = n+1, data = val$x, byrow = F)
    zmat <- matrix(ncol = n + 1, nrow = n + 1, pcopula(object, xmat, ymat), byrow = F)
    val$z <- zmat
#    contour(val,nlevels=nlevels,add=add, ..., xlab = "u",ylab= "v")
    contour(x=c(0:n)/n,y=c(0:n)/n,z=val$z, nlevels = nlevels, add = add, ..., xlab = "u", ylab = "v")
    invisible(val)
} )


#############################################################
# persp functions for 3-D plots

## persp.dbivd <- function(dist, n = 50, xlim = NA, ylim = NA, ...) {

##     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
##     divis <- divis[2:(n+1)]

##     if ( (is.na(xlim[1])) | (length(xlim) != 2)) {  
##       x <- evalFunc(divis,get(paste("q",dist@Xmarg,sep = "")),dist@param.Xmarg)
##     } else {  x <- seq(from = xlim[1], to = xlim[2], length = n) }

##     if (is.na(ylim[1]) | length(ylim) != 2) {
##      y <- evalFunc(divis,get(paste("q",dist@Ymarg,sep = "")),dist@param.Ymarg)
##     } else { y <- seq(from = ylim[1], to = ylim[2], length = n) }

##     xmat <- rep(x, n )
##     ymat <- rep(y, each = n )
##     zmat <- dbivd(dist, xmat, ymat)

## persp3d(x, x, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "Density")

## #    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat),
## #             z = as.vector(zmat)))
## #    guiPlot("Data Grid Surface", DataSetValues  = val2, ...) 

##     val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, 
##        data = zmat   ))
##   invisible(val)     

## }

## persp.pbivd <- function(dist, n = 50, xlim = NA, ylim = NA, ...) {

##     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
##     divis <- divis[2:(n+1)]

##     if ( (is.na(xlim[1])) | (length(xlim) != 2)) {  
##       x <- evalFunc(divis,get(paste("q",dist@Xmarg,sep = "")),dist@param.Xmarg)
##     } else {  x <- seq(from = xlim[1], to = xlim[2], length = n) }

##     if (is.na(ylim[1]) | length(ylim) != 2) {
##      y <- evalFunc(divis,get(paste("q",dist@Ymarg,sep = "")),dist@param.Ymarg)
##     } else { y <- seq(from = ylim[1], to = ylim[2], length = n) }

##     xmat <- rep(x, n )
##     ymat <- rep(y, each = n )
##     zmat <- pbivd(dist, xmat, ymat)

## persp3d(x, y, zmat, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "cdf")

## #    val2 <- data.sheet(list(x = as.vector(xmat), y = as.vector(ymat),
## #             z = as.vector(zmat)))   
## #    guiPlot("Data Grid Surface", DataSetValues  = val2, ...) 
##  val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, 
##        data = zmat   ))
##   invisible(val)     




## }

## contour.dbivd <- function(dist, n = 50, xlim = NA, ylim = NA, ...) {

##     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
##     divis <- divis[2:(n+1)]

##     if ( (is.na(xlim[1])) | (length(xlim) != 2)) {  
##       x <- evalFunc(divis,get(paste("q",dist@Xmarg,sep = "")),dist@param.Xmarg)
##     } else {  x <- seq(from = xlim[1], to = xlim[2], length = n) }

##     if (is.na(ylim[1]) | length(ylim) != 2) {
##      y <- evalFunc(divis,get(paste("q",dist@Ymarg,sep = "")),dist@param.Ymarg)
##     } else { y <- seq(from = ylim[1], to = ylim[2], length = n) }

##     xmat <- rep(x, n )
##     ymat <- rep(y, each = n )
##     zmat <- dbivd(dist, xmat, ymat)
##   val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, 
##        data = zmat   ))
## contour(x,y,val$z)
## title("Contour Plot of the Density", xlab="u",ylab="v")
##   invisible(val)


## }

## contour.pbivd <- function(dist, n = 50, xlim = NA, ylim = NA, ...) {

##     divis <- seq(from = 0.001, to = 0.999, length = (n+2)) 
##     divis <- divis[2:(n+1)]

##     if ( (is.na(xlim[1])) | (length(xlim) != 2)) {  
##       x <- evalFunc(divis,get(paste("q",dist@Xmarg,sep = "")),dist@param.Xmarg)
##     } else {  x <- seq(from = xlim[1], to = xlim[2], length = n) }

##     if (is.na(ylim[1]) | length(ylim) != 2) {
##      y <- evalFunc(divis,get(paste("q",dist@Ymarg,sep = "")),dist@param.Ymarg)
##     } else { y <- seq(from = ylim[1], to = ylim[2], length = n) }

##     xmat <- rep(x, n )
##     ymat <- rep(y, each = n )
##     zmat <- pbivd(dist, xmat, ymat)

##    val <- list(x = x, y = y, z = matrix(ncol = n, nrow = n, byrow = F, 
##        data = zmat   ))
## contour(x,y,val$z)
## title("Contour Plot of the CDF", xlab="u",ylab="v")
##      invisible(val)
## }


######################### 
### Functions lambda  ### 

setMethod("lambda","frank.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","kimeldorf.sampson.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","bb1.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","bb2.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","bb3.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","bb6.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","bb7.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )

setMethod("lambda","joe.copula",
     function(copula, t) {
    Phi(copula, t)/PhiDer(copula, t)
} )



##################################### 
### Method Functions dcdx
#####################################


setMethod("dcdx","gumbel.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    lx <- log(x)
    ly <- log(y)
    lxy <- lx + ly
    t <- lx/lxy
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    val[good] <- 1 - (exp(ad0 * lxy) * (ad0 * lxy + ad1 * ( - lx + lxy)))/(x * lxy)
    val[u == 0 | u == 1] <- 1
    val[v == 0 & (u >= 0) & (u <= 1)] <- 1
    val[v == 1 & (u >= 0) & (u <= 1)] <- 0
    val
} )

setMethod("dcdx","galambos.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    lx <- log(x)
    ly <- log(y)
    lxy <- lx + ly
    t <- lx/lxy
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    val[good] <- 1 - (exp(ad0 * lxy) * (ad0 * lxy + ad1 * ( - lx + lxy)))/(x * lxy)
    val[u == 0 | u == 1] <- 1
    val[v == 0 & (u >= 0) & (u <= 1)] <- 1
    val[v == 1 & (u >= 0) & (u <= 1)] <- 0
    val
} )

setMethod("dcdx","husler.reiss.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    lx <- log(x)
    ly <- log(y)
    lxy <- lx + ly
    t <- lx/lxy
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    val[good] <- 1 - (exp(ad0 * lxy) * (ad0 * lxy + ad1 * ( - lx + lxy)))/(x * lxy)
    val[u == 0 | u == 1] <- 1
    val[v == 0 & (u >= 0) & (u <= 1)] <- 1
    val[v == 1 & (u >= 0) & (u <= 1)] <- 0
    val
} )

setMethod("dcdx","tawn.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    lx <- log(x)
    ly <- log(y)
    lxy <- lx + ly
    t <- lx/lxy
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    val[good] <- 1 - (exp(ad0 * lxy) * (ad0 * lxy + ad1 * ( - lx + lxy)))/(x * lxy)
    val[u == 0 | u == 1] <- 1
    val[v == 0 & (u >= 0) & (u <= 1)] <- 1
    val[v == 1 & (u >= 0) & (u <= 1)] <- 0
    val
} )

setMethod("dcdx","bb5.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    lx <- log(x)
    ly <- log(y)
    lxy <- lx + ly
    t <- lx/lxy
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    val[good] <- 1 - (exp(ad0 * lxy) * (ad0 * lxy + ad1 * ( - lx + lxy)))/(x * lxy)
    val[u == 0 | u == 1] <- 1
    val[v == 0 & (u >= 0) & (u <= 1)] <- 1
    val[v == 1 & (u >= 0) & (u <= 1)] <- 0
    val
} )


setMethod("dcdx","frank.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","kimeldorf.sampson.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","bb1.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","bb2.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","bb3.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","bb6.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","bb7.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

setMethod("dcdx","joe.copula",
     function(copula,u,v) {
    val <- vector(length = length(u))
    good <- u > 0 & u < 1
    x <- u[good]
    y <- v[good]
    iphip <- InvPhifirstDer(copula, Phi(copula, x) + Phi(copula, y))
    phip <- PhiDer(copula, x)
    val[good] <- iphip * phip
    val[u == 0] <- 1
    val[u == 1] <- v
    val
} )

dcdx.ev.copula <- function(copula,u,v) {
   val <- vector (length = length(u))
    good <- u > 0 & u < 1 & v > 0 & v < 1
    x <- u[good]
    y <- v[good]
    lx <- log(x)
    ly <- log(y)
    lxy <- lx +ly

    t <- lx/lxy
    ad0 <- A(copula, t)
    ad1 <- AfirstDer(copula, t)
    val[good] <- 1 - (exp(ad0*lxy)*(ad0*lxy + ad1*(-lx + lxy)))/(x*lxy)

   val[u == 0 | u == 1] <- 1    
   val[v == 0 & (u >= 0) & (u <= 1)] <- 1   
   val[v == 1 & (u >= 0) & (u <= 1)] <- 0        
   val
}

dcdx.archm.copula <- function(copula,u,v) {
  val <- vector (length = length(u))
    good <- u > 0 & u < 1

    x <- u[good]
    y <- v[good]

    iphip <- InvPhifirstDer(copula, Phi(copula,x)+Phi(copula,y))
    phip <- PhiDer(copula,x)
    val[good] <- iphip * phip   

   val[u == 0] <- 1
   val[u == 1] <- v     
   val
}


################################
# VaR Examples Computations

fbar.exp.portf <- function(a, copula, x.est, y.est, lambda1, lambda2, q = NA,
       subdivisions = 1000, tol.f = 5e-5) 
{
    a <- -a
    BE <- exp(a)/lambda2
    UL <- gpd.2p(a - log(lambda1) , x.est)
    ftem <- function(x,BE,copula, x.est, y.est, lambda1, lambda2)
    {
        val <- vector(length = length(x))
        u <- gpd.2q(x,x.est)
        uuu <- BE - lambda1/lambda2*exp(u)
        val[uuu <= 0] <- 1
        val[uuu >= exp(1)] <- 1
        v <- log(uuu[uuu>0 & uuu <exp(1)])
        px <- gpd.2p(v , y.est)
       val[uuu>0 & uuu <exp(1) ] <- dcdx(copula,x[uuu>0 & uuu <exp(1)] ,px)
     val
    }

   ret <- integrate(ftem, 0, UL, abs.tol= tol.f, subdivisions = subdivisions, 
BE=BE,copula=copula,x.est=x.est,y.est=y.est,lambda1=lambda1,lambda2=lambda2)$value
   if (!is.na(q)) ret <- ret - q
   ret
}

fbar.lin.portf <- function(Q, copula, x.est, y.est, lambda1, lambda2, 
        subdivisions = 1000, tol.f = 5e-5)  
{   
    ftem <- function(x,Q,copula, x.est, y.est, lambda1, lambda2)
    {
        qq <- (Q - gpd.1q(x,y.est)*lambda2)/lambda1
        px <- vector(length = length(qq))
        px[is.inf(qq) & qq < -1] <- 0
        px[is.inf(qq) & qq > 1] <- 1
        px[!is.inf(qq) ] <- gpd.1p(qq[!is.inf(qq) ] , x.est)
        val <- vector(length = length(x))
        val[ px >=1] <- 1
        val[ px <= 0] <- 1
        val[px > 0 & px < 1] <- dcdx.max.copula(copula,x[px > 0 & px < 1] ,px[px > 0 & px < 1] )
        val
    }
   integrate(ftem, 0, 1,Q = Q,copula = copula, x.est = x.est, y.est = y.est,
       lambda1 = lambda1, lambda2 = lambda2, subdivisions = 1000, abs.tol = 5e-5)$value
}

VaR.exp.portf <- function(Q, copula, x.est, y.est, lambda1, lambda2,
      range, subdivisions = 1000, tol.f = 5e-5)   
{     
     uniroot(fbar.exp.portf, interval=range, copula = copula,x.est = x.est, y.est = y.est,
       lambda1 = lambda1, lambda2 = lambda2, q = Q,
          subdivisions = subdivisions,tol.f=tol.f)$root
}

VaR.exp.sim <- function(n, Q, copula, x.est, y.est, lambda1, lambda2 ) 
{
    sim.uv <- rcopula(copula, n)
    Xsim <- gpd.2q(sim.uv$x,x.est)
    Ysim <- gpd.2q(sim.uv$y,y.est)
    ww <- log(lambda1*exp(Xsim) + lambda2*exp(Ysim))
    varv <- -as.vector(quantile(ww,Q))
    esf <- vector(length = length(Q), mode = "numeric")
    for (i in 1:length(Q))  esf[i] <- -mean(ww[ww < -varv[i]])

    val <- c(n,varv,esf)
    nn <-  vector(length = length(Q)*2 + 1, mode = "character")
    nn[1] <- "Simulation size"
    nn[2:(length(Q)+1)] <- paste("VaR Q=",Q,sep = "")
    nn[(length(Q)+2):(2*length(Q)+1)] <- paste("ES Q=",Q,sep = "")
    names(val) <- nn
    val
}


#########################
# Kendall Table

"Kend.table" <- 
structure(.Data = list("tau" = c(0., 0.10000000000000001, 0.20000000000000001, 0.29999999999999999, 0.40000000000000002,
    0.5, 0.59999999999999998, 0.69999999999999996, 0.80000000000000004, 
    0.90000000000000002, 1.)
, "B1" = c(0., 0.1564344650402309, 0.3090169943749474, 0.45399049973954669, 0.58778525229247314,
    0.70710678118654746, 0.80901699437494745, 0.89100652418836779, 0.95105651629515364,
    0.98768834059513777, 1.)
, "B3" = c(0., 0.91000000000000003, 1.8600000000000001, 2.9199999999999999, 4.1600000000000001,
    5.7400000000000002, 7.9299999999999997, 11.4, 18.199999999999999, 
    20.899999999999999, Inf)
, "B4" = c(0., 0.22, 0.5, 0.85999999999999999, 1.3300000000000001, 2., 3., 4.6699999999999999, 8.,
    18., Inf)
, "B6" = c(1., 1.1111111111111109, 1.25, 1.4285714285714279, 1.666666666666667, 2., 2.5, 
    3.333333333333333, 5.0000000000000009, 10., Inf)
, "B7" = c(0., 0.34000000000000002, 0.51000000000000001, 0.69999999999999996, 0.94999999999999996,
    1.28, 1.79, 2.6200000000000001, 4.29, 9.3000000000000007, Inf)
, "B8" = c(0., 0.66000000000000003, 0.87, 1.1100000000000001, 1.4099999999999999, 1.8, 
    2.3900000000000001, 3.3399999999999999, 5.2400000000000002, 10.9, Inf)
, "BB1.delta" = c(1., 1.1000000000000001, 1.2, 1.3, 1.3999999999999999, 1.6499999999999999, 2., 
    2.6499999999999999, 3., 5., Inf)
, "bb1.theta" = c(0., 0.020202020202019891, 0.083333333333333481, 0.1978021978021976, 0.38095238095238099,
    0.42424242424242431, 0.5, 0.51572327044025146, 1.3333333333333339, 
    2.0000000000000009, Inf)
, "Twan.copula" = c(0., 1.1258349999999999, 1.290394, 1.5154939999999999, 1.843656, 2.3701919999999999, 
    3.3637899999999998, 5.997922, 37.203530000000001, 172.10230000000001, Inf)
, "TC.a.b" = c(0.90000000000000002, 0.90000000000000002, 0.90000000000000002, 0.90000000000000002, 
    0.90000000000000002, 0.90000000000000002, 0.90000000000000002, 0.90000000000000002,
    0.90000000000000002, 0.94999999999999996, Inf)
, "BB2.theta" = c(0., 0.14999999999999999, 0.20000000000000001, 0.29999999999999999, 0.40000000000000002,
    0.5, 0.59999999999999998, 0.69999999999999996, 0.78000000000000003, 
    0.81000000000000005, Inf)
, "BB2.delta" = c(0., 0.41547420000000002, 1.2451000000000001, 1.4251780000000001, 1.6769670000000001,
    2.0480670000000001, 2.6324589999999999, 3.6496490000000001, 5.9888640000000004,
    8.5999999999999996, Inf)
, "BB3.theta" = c(1., 1.05, 1.1000000000000001, 1.2, 1.3, 1.3999999999999999, 1.5, 1.8, 2.1000000000000001,
    2.2010000000000001, Inf)
, "BB3.delta" = c(0., 0.1152268, 0.26838640000000002, 0.37031439999999999, 0.54731560000000001, 
    0.84593099999999999, 1.37849, 1.8886909999999999, 3.7887529999999998, 
    13.029999999999999, Inf)
, "BB5.theta" = c(1., 1.05, 1.1000000000000001, 1.2, 1.3, 1.3999999999999999, 1.5, 1.8, 1.8999999999999999,
    1.95, Inf)
, "BB5.delta" = c(0., 0.26372449999999997, 0.37258530000000001, 0.43810700000000002, 0.54184339999999998,
    0.69927689999999998, 0.94607929999999996, 1.1348119999999999, 1.9201589999999999,
    4.4188349999999996, Inf)
, "BB6.theta" = c(1., 1.05, 1.1000000000000001, 1.1499999999999999, 1.2, 1.25, 1.3500000000000001, 1.55,
    1.8, 2., Inf)
, "BB6.delta" = c(1., 1.079931, 1.18201, 1.31548, 1.495727, 1.750732, 2.0865680000000002, 
    2.5480239999999998, 3.4635950000000002, 6.4493419999999997, Inf)
, "BB7.theta" = c(1., 1.05, 1.1000000000000001, 1.25, 1.5, 1.8, 2.2999999999999998, 2.7999999999999998,
    3.5, 9.5, Inf)
, "BB7.delta" = c(0., 0.1649669, 0.3882816, 0.58719399999999999, 0.81623789999999996, 1.2346159999999999,
    1.9072309999999999, 3.8238219999999998, 10.30184, 60.21819, Inf)
, "BB4.theta" = c(0., 0.19, 0.29999999999999999, 0.40000000000000002, 0.5, 0.69999999999999996, 
    0.90000000000000002, 1.3, 1.8999999999999999, 2.2999999999999998, Inf)
, "BB4.delta" = c(0., 0.1750353, 0.30702400000000002, 0.43810700000000002, 0.59777539999999996, 
    0.75476860000000001, 1.004845, 1.3052299999999999, 1.8524259999999999, 
    3.9417589999999998, Inf)
)
, names = c("tau", "B1", "B3", "B4", "B6", "B7", "B8", "BB1.delta", "bb1.theta", "Twan.copula", 
    "TC.a.b", "BB2.theta", "BB2.delta", "BB3.theta", "BB3.delta", "BB5.theta", 
    "BB5.delta", "BB6.theta", "BB6.delta", "BB7.theta", "BB7.delta", "BB4.theta", 
    "BB4.delta")
, row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11")
, .guiColInfo = list(c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-2")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-3")
, c("TXPROP_ColWidth<-12", "TXPROP_ColJustification<-Right", "TXPROP_ObjectNote<-", 
    "TXPROP_ColFloatFormat<-Decimal", "TXPROP_ColPrecision<-5")
)
, class = "data.frame"
)


#################
# bivd functions



"bivd" <-  function(cop, marginX = "unif", marginY = marginX, param.marginX = c(0, 1),
    param.marginY = param.marginX)
{
    yyy <- class(cop)
    cl <- substring(yyy, 1, nchar(yyy) - 7)
    val <- new(paste(cl, "bivd", sep = "."), copula = cop, Xmarg = marginX,
        Ymarg = marginY, param.Xmarg = param.marginX, param.Ymarg =
        param.marginY)
    val
}

## "persp.dbivd" <-  function(dist, n = 50, xlim = NA, ylim = NA, ...)
## {
##     divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
##     divis <- divis[2:(n + 1)]
##     if((is.na(xlim[1])) | (length(xlim) != 2)) {
##         x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")),
##             dist@param.Xmarg)
##     }
##     else {
##         x <- seq(from = xlim[1], to = xlim[2], length = n)
##     }
##     if(is.na(ylim[1]) | length(ylim) != 2) {
##         y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")),
##             dist@param.Ymarg)
##     }
##     else {
##         y <- seq(from = ylim[1], to = ylim[2], length = n)
##     }
##     xmat <- rep(x, n)
##     ymat <- rep(y, each = n)
##     zmat <- dbivd(dist, xmat, ymat)
## #    print(wireframe(zmat ~ xmat * ymat, aspect = c(1, 8.5/11), scales =
## #        list(arrows = F), zlab = "Density", zoom = 1.2, xlab = "x",
## #        ylab = "y"))
##    val <- list(x = x, y = y, z = matrix(zmat, n, n))
## persp3d(x, y, val$z, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "Density")
##     invisible(val)
## }

## "persp.pbivd" <- function(dist, n = 50, xlim = NA, ylim = NA, ...)
## {
##     divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
##     divis <- divis[2:(n + 1)]
##     if((is.na(xlim[1])) | (length(xlim) != 2)) {
##         x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")), dist@
##             param.Xmarg)
##     }
##     else {
##         x <- seq(from = xlim[1], to = xlim[2], length = n)
##     }
##     if(is.na(ylim[1]) | length(ylim) != 2) {
##         y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")), dist@
##             param.Ymarg)
##     }
##     else {
##         y <- seq(from = ylim[1], to = ylim[2], length = n)
##     }
##     xmat <- rep(x, n)
##     ymat <- rep(y, each = n)
##     zmat <- pbivd(dist, xmat, ymat)
##     val <- list(x = x, y = y, z = matrix(zmat, n, n))
## persp3d(x, y, val$z, aspect=c(1, 1, 0.5), col = "lightblue",xlab = "x", ylab = "y", zlab = "cdf")

## #    print(wireframe(zmat ~ xmat * ymat, aspect = c(1, 8.5/11), scales = list(arrows
## #         = F), zlab = "Density", zoom = 1.2, xlab = "x", ylab = "y"))
##     invisible(val)
## }


## "contour.dbivd" <- function(dist, n = 50, xlim = NA, ylim = NA, ...)
## {
##     divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
##     divis <- divis[2:(n + 1)]
##     if((is.na(xlim[1])) | (length(xlim) != 2)) {
##         x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")),
##             dist@param.Xmarg)
##     }
##     else {
##         x <- seq(from = xlim[1], to = xlim[2], length = n)
##     }
##     if(is.na(ylim[1]) | length(ylim) != 2) {
##         y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")),
##             dist@param.Ymarg)
##     }
##     else {
##         y <- seq(from = ylim[1], to = ylim[2], length = n)
##     }
##     xmat <- rep(x, n)
##     ymat <- rep(y, each = n)
##     zmat <- dbivd(dist, xmat, ymat)
##     print(contourplot(zmat ~ xmat * ymat, xlab = "x", ylab = "y", cuts =
##         10))
##     val <- list(x = x, y = y, z = matrix(zmat, n, n))
##     invisible(val)
## }

## "contour.pbivd" <- function(dist, n = 50, xlim = NA, ylim = NA, ...)
## {
##     divis <- seq(from = 0.001, to = 0.999, length = (n + 2))
##     divis <- divis[2:(n + 1)]
##     if((is.na(xlim[1])) | (length(xlim) != 2)) {
##         x <- evalFunc(divis, get(paste("q", dist@Xmarg, sep = "")),
##             dist@param.Xmarg)
##     }
##     else {
##         x <- seq(from = xlim[1], to = xlim[2], length = n)
##     }
##     if(is.na(ylim[1]) | length(ylim) != 2) {
##         y <- evalFunc(divis, get(paste("q", dist@Ymarg, sep = "")),
##             dist@param.Ymarg)
##     }
##     else {
##         y <- seq(from = ylim[1], to = ylim[2], length = n)
##     }
##     xmat <- rep(x, n)
##     ymat <- rep(y, each = n)
##     zmat <- pbivd(dist, xmat, ymat)
##     val <- list(x = x, y = y, z = matrix(zmat, n, n))
##     print(contourplot(zmat ~ xmat * ymat, xlab = "x", ylab = "y", cuts =
##         10))
##     invisible(val)
## }



"gpdjoint.1p" <- function(x, y, copula, x.est, y.est)
{
    h <- gpd.1p(x, x.est)
    g <- gpd.1p(y, y.est)
    val <- pcopula(copula, h, g)
    val
}

"gpdjoint.2p" <- function(x, y, copula, x.est, y.est)
{
    h <- gpd.2p(x, x.est)
    g <- gpd.2p(y, y.est)
    val <- pcopula(copula, h, g)
    val
}


##################
# Stand alone (default) empirical.copula function

empirical.copula <- function(x, y = NA) 
{
    data <- NA
    if(inherits(x,"empirical.copula")) {data <- x}
    else if ( is.list(x))
    {
       temp <- data.frame(x)
       data <- new("empirical.copula",x = temp[,1], y = temp[,2])
    }
    else if (!is.na(y[1]) & length(x) == length(y))
    {
       data <- new("empirical.copula",x = x, y = y)
    }
    else
    {
      stop("Can't create empirical copula. Unknown parameters supplied")
    }

    inputsOK <- T
    message <- ""
    if (sum(data@x < 0) > 0 | sum(data@y < 0) > 0) 
    {
        inputsOK <- F
        message <-
           paste(message," At least one of the observed data points has a negative coordinate", sep = " \n")
    }
    if (sum(data@x > 1) > 0 | sum(data@y > 1) > 0) 
    {
        inputsOK <- F
        message <-
           paste(message," At least one of the observed data points has a coordinate greater than 1",
                  sep = " \n")
    }
    if (!inputsOK) {stop(message)}         
    data    
}

Try the Rsafd package in your browser

Any scripts or data that you put into this service are public.

Rsafd documentation built on May 2, 2019, 5:20 p.m.