# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General 
# Public License along with this library; if not, write to the 
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
# MA  02111-1307  USA

# Copyrights (C)
# for this R-port: 
#   1999 - 2007, Diethelm Wuertz, GPL
#   Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
#   info@rmetrics.org
#   www.rmetrics.org
# for the code accessed (or partly included) from other R-ports:
#   see R's copyright and license files
# for the code accessed (or partly included) from contributed R-ports
# and other sources
#   see Rmetrics's copyright file


################################################################################
# FUNCTION:                  ELLIPTICAL COPULAE RANDOM DEVIATES:
#  rellipticalCopula          Generates elliptical copula variates
#  rellipticalSlider          Generates interactive plots of random variates
#  .rnormCopula               Generates normal copula random variate
#  .rcauchyCopula             Generates Cauchy copula random variate
#  .rtCopula                  Generates Student-t copula random variate
# FUNCTION:                  ELLIPTICAL COPULAE PROBABILITY:
#  pellipticalCopula          Computes elliptical copula probability
#  pellipticalSlider          Generates interactive plots of probability
#  .pnormCopula               Computes normal copula probability
#  .pcauchyCopula             Computes Cauchy copula probability
#  .ptCopula                  Computes Student-t copula probability
#  .pellipticalCopulaGrid     Fast equidistant grid version
#  .pellipticalCopulaDiag     Fast diagonal cross section version
#  .pellipticalPerspSlider    Interactive perspective plots of probability
#  .pellipticalContourSlider  Interactive contour plots of probability
# FUNCTION:                  ELLIPTICAL COPULAE DENSITY:
#  dellipticalCopula          Computes elliptical copula density 
#  dellipticalSlider          Generates interactive plots of density
#  .dnormCopula               Computes normal copula density
#  .dcauchyCopula             Computes Cauchy copula density
#  .dtCopula                  Computes Student-t copula density
#  .dellipticalCopulaGrid     Fast grid version for elliptical copula density
#  .dellipticalPerspSlider    Interactive perspective plots of density
#  .dellipticalContourSlider  Interactive contour plots of density
################################################################################


################################################################################
# FUNCTION:                  ELLIPTICAL COPULAE RANDOM DEVIATES:
#  rellipticalCopula          Generates elliptical copula variates
#  rellipticalSlider          Generates interactive plots of random variates
#  .rnormCopula               Generates normal copula random variate
#  .rcauchyCopula             Generates Cauchy copula random variate
#  .rtCopula                  Generates Student-t copula random variate


rellipticalCopula =
function(n, rho = 0.75, param = NULL, type = c("norm", "cauchy", "t"))
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes extreme value copula probability
    
    # Arguments:
    #   n - number of deviates to be generated.
    #   rho - a numeric value setting the coorelation strength, ranging
    #       between minus one and one.
    #   nu - the number of degrees of freedom, only required for 
    #       Student-t copulae.
    #   type - the type of the elliptical copula. Either "norm" or
    #       "t" denoting the normal or Student-t copula, respectively.
    #   output - a character string specifying how the output should
    #       be formatted. By default a vector of the same length as 
    #       'u' and 'v'. If specified as "list" then 'u' and 'v' are
    #       expected to span a two-dimensional grid as outputted by the 
    #       function 'grid2d' and the function returns a list with
    #       elements '$x', 'y', and 'z' which can be directly used 
    #       for example by 2D plotting functions. 
    
    # Value:
    #   returns a vector or list of probabilities depending on the
    #   value of the "output" variable.
    
    # Example:
    #   Diagonal Value: pnormCopula((0:10)/10)
    #   persp(pnormCopula(u = grid2d(), output = "list"))
    
    # FUNCTION:
    
    # Settings:
    type = match.arg(type)
    
    # Parameters:
    if (type == "t") {
        if(is.null(param)) {
            param = c(nu = 4)
        } else {
            param = c(nu = param)
        }
        names(param) = "nu"
    }
    
    # Copula:
    if (type == "norm") 
        ans = .rnormCopula(n = n, rho = rho)
    if (type == "cauchy") 
        ans = .rcauchyCopula(n = n, rho = rho)
    if (type == "t") 
        ans = .rtCopula(n = n, rho = rho, nu = param)
    
    # Add Control Attribute:
    control = list(rho = rho, param = param, type = type)
    attr(ans, "control")<-unlist(control)
        
    # Return Value:
    ans
}


# ------------------------------------------------------------------------------


rellipticalSlider = 
function(B = 100)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively perspective plots of random variates
    
    #FUNCTION:
    
    # Graphic Frame:
    par(mfrow = c(1, 1))
    
    # Internal Function:
    refresh.code = function(...)
    {
        # Startup Counter:
        .counter <<- .counter + 1
        if (.counter < 7) return ()
        
        # Sliders:
        Copula = .sliderMenu(no = 1)
        N = .sliderMenu(no = 2)
        rho = .sliderMenu(no = 3)
        nu = .sliderMenu(no = 4)
        seed = .sliderMenu(no = 5)
        size = .sliderMenu(no = 6)
        col = .sliderMenu(no = 7)
        Names = c("- Normal", "- Cauchy", "- Student t") 
        Type = c("norm", "cauchy", "t")
        eps = 1.0e-6
        if (rho == +1) rho = rho - eps
        if (rho == -1) rho = rho + eps
        
        # Tau and Rho:
        Tau = ellipticalTau(rho)
        Rho = ellipticalRho(rho)
        
        # Plot: 
        Title = paste("Elliptical Copula No:", as.character(Copula), 
            Names[Copula], "\nrho =", as.character(rho), "|")
        if (Copula == 2) Title = paste(Title, "nu =", as.character(nu), "|")   
        Title = paste(Title,       
            "Kendall = ", as.character(round(Tau, digits = 3)), "|",
            "Spearman = ", as.character(round(Rho, digits = 3)) )       
        set.seed(seed)
        R = rellipticalCopula(n = N, rho = rho, param = nu, type = Type[Copula])
        plot(x = R[, 1], y = R[, 2], xlim = c(0, 1), ylim = c(0, 1), 
            xlab = "u", ylab = "v", pch = 19, col = col, cex = size)
        title(main = Title)
        
        # Reset Frame:
        par(mfrow = c(1, 1))
    }
  
    # Open Slider Menu:
    .counter <<- 0 
    plot.names = c("Plot - size", "... color")
    .sliderMenu(refresh.code,
        names       = c("Copula",   "N", "rho", "t: nu", "seed", plot.names),
        minima      = c(       1,  1000,    -1,       1,   1000,     0,   1),
        maxima      = c(       3, 10000,    +1,       B,   9999,     1,  16),
        resolutions = c(       1,   500,  0.01,       1,      1,   0.1,   1),
        starts      = c(       1,  1000,     0,       4,   4711,   0.5,   1)) 
}


# ------------------------------------------------------------------------------


.rnormCopula =
function(n, rho = 0.75)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Generates normal copula random variate
    
    # Example:
    #   UV = rnormCopula(n = 10000); plot(UV[,1], UV[,2], cex = 0.25)

    # FUNCTION:
    
    # Use: X = .rnorm2d(n, rho) or alternatively:
    X = .rnorm2d(n = n, rho = rho)
    
    # Generate
    Z <- NULL
    for(i in (1:n)) Z <- rbind(Z, pnorm(X [i,]))
    
    # Return Value:
    Z
}


# ------------------------------------------------------------------------------


.rcauchyCopula =
function(n, rho = 0.75)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Generates Student-t copula random variate
    
    # Example:
    #   UV = rtCopula(n = 10000); plot(UV[,1], UV[,2], cex = 0.25)

    # FUNCTION:
    
    # Cauchy Deviates:
    Z = .rtCopula(n = n, rho = rho, nu = 1)
    
    # Return Value:
    Z
}


# ------------------------------------------------------------------------------


.rtCopula =
function(n, rho = 0.75, nu = 4)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Generates Student-t copula random variate
    
    # Example:
    #   UV = rtCopula(n = 10000); plot(UV[,1], UV[,2], cex = 0.25)

    # FUNCTION:
    
    # Use: X = .rnorm2d(n, rho) or alternatively:
    X = rt2d(n = n, rho = rho, nu = nu)
    
    # Generate
    Z = NULL
    for (i in (1:n)) Z = rbind(Z, pt(X [i,], df = nu))
    
    # Return Value:
    Z
}


################################################################################
# FUNCTION:                  ELLIPTICAL COPULAE PROBABILITY:
#  pellipticalCopula          Computes elliptical copula probability
#  pellipticalSlider          Generates interactive plots of probability
#  .pnormCopula               Computes normal copula probability
#  .pcauchyCopula             Computes Cauchy copula probability
#  .ptCopula                  Computes Student-t copula probability
#  .pellipticalCopulaGrid     Fast equidistant grid version
#  .pellipticalCopulaDiag     Fast diagonal cross section version
#  .pellipticalPerspSlider    Interactive perspective plots of probability
#  .pellipticalContourSlider  Interactive contour plots of probability


pellipticalCopula =
function(u = 0.5, v = u, rho = 0.75, param = NULL, type = ellipticalList(), 
output = c("vector", "list"), border = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes extreme value copula probability
    
    # Arguments:
    #   u, v - two numeric values or vectors of the same length at
    #       which the copula will be computed. If 'u' is a list then the
    #       the '$x' and '$y' elements will be used as 'u' and 'v'.
    #       If 'u' is a two column matrix then the first column will
    #       be used as 'u' and the the second as 'v'.
    #   rho - a numeric value setting the coorelation strength, ranging
    #       between minus one and one.
    #   param - distributional parameters, the number of degrees of 
    #       freedom for the Student-t copulae.
    #   type - the type of the elliptical copula. Either "norm" or
    #       "t" denoting the normal or Student-t copula, respectively.
    #   output - a character string specifying how the output should
    #       be formatted. By default a vector of the same length as 
    #       'u' and 'v'. If specified as "list" then 'u' and 'v' are
    #       expected to span a two-dimensional grid as outputted by the 
    #       function 'grid2d' and the function returns a list with
    #       elements '$x', 'y', and 'z' which can be directly used 
    #       for example by 2D plotting functions. 
    
    # Value:
    #   returns a vector or list of probabilities depending on the
    #   value of the "output" variable.
    
    # FUNCTION:

    # Match Arguments:
    type = match.arg(type)
    output = match.arg(output)
    
    # Settings:
    subdivisions = 100
    if (is.list(u)) {
        v = u[[2]]
        u = u[[1]]
    }
    if (is.matrix(u)) {
        v = u[, 2]
        u = u[, 1]
    }
    if (length(u) == 1 & u[1] > 1) {
        return(.pellipticalCopulaGrid(N = u, rho, param, type, border = border))
    }
    
    # Parameters:
    if (type == "t") if (is.null(param)) param = c(nu = 4)
    if (type == "kotz") if (is.null(param)) param = c(r = 1)
    if (type == "epower") if (is.null(param)) param = c(r = 1, s = 1)
    
    # Specical Copulae:
    if (type == "norm") {
        if (rho == -1) {
            ans = pfrechetCopula(u = u, v = v, type = "m", output = output)
            return(ans)
        } else if (rho == +1) {
            ans = pfrechetCopula(u = u, v = v, type = "w", output = output)
            return(ans)
        } else {
            ans = .pnormCopula(u = u, v = v, rho = rho, output = output)
            return(ans)
        }
    } else if (type == "cauchy") {
        ans = .pcauchyCopula(u = u, v = v, rho = rho, output = output)
        return(ans)
    } else if (type == "t") {
        if (is.null(param)) param = 4
        ans = .ptCopula(u = u, v = v, rho = rho, nu = param, output = output)
        return(ans)
    }
    
    # The remaining Copulae - Compute Density on Regular Grid:
    N = subdivisions
    x = (0:N)/N
    c.uv = .dellipticalCopulaGrid(N = N, rho, param, type, border = TRUE)
    c.uv$z[is.na(c.uv$z)] = 0      
    
    # Integrate to get Probability:
    C.uv = 0*c.uv$z
    for (i in 1:(N+1)) {
        D = matrix(rep(0, times = (N+1)^2), ncol = N+1)
        for (j in 1:i) {    
            D[1:i, j] = 1
            C.uv[i,j] = C.uv[j,i] = sum(D*c.uv$z)
        }
    }
    C.uv = C.uv/N^2
    
    # Take care about the Boundary on the Unit Square:
    C.uv[1, ] = C.uv[, 1] = 0
    C.uv[N+1, ] = C.uv[, N+1] = c.uv$x   
    
    # Interpolate for the desired Values on the grid:
    U0 = trunc(u*N)
    V0 = trunc(v*N)
    P = (u - U0/N)
    Q = (v - V0/N)
    U0 = U0 + 1
    U1 = U0 + 1
    V0 = V0 + 1
    V1 = V0 + 1
    C.vec = rep(NA, times = length(u))
    for ( i in 1:length(u) ) {
        p = P[i]
        q = Q[i]
        if (p == 0 & q == 0) {
            C.vec[i] = C.uv[U0[i], V0[i]] 
        } else if (p == 0 & q > 0) {
            C.vec[i] = (1-q)*C.uv[U0[i], V0[i]] + q*C.uv[U0[i], V1[i]]
        } else if (p > 0 & q == 0) {
            C.vec[i] = (1-p)*C.uv[U0[i], V0[i]] + p*C.uv[U1[i], V0[i]] 
        } else {
            C.vec[i] = (1-p)*(1-q)*C.uv[U0[i], V0[i]] + 
                p*(1-q)*C.uv[U1[i], V0[i]] + (1-p)*q*C.uv[U0[i], V1[i]] + 
                p*q*C.uv[U1[i], V1[i]] 
        }
    }
    C.uv = round(C.vec, digits = 3)
    attr(C.uv, "control") <- c(rho = rho)
    
    # As List ?
    if (output == "list") {
        N = sqrt(length(u))
        x = u[1:N]
        y = matrix(v, ncol = N)[1, ]
        names(x) = names(y) = NULL
        C.uv = list(x = x, y = y,  z = matrix(C.uv, ncol = N))
    }
    
    # Return Value:
    C.uv 
}


# ------------------------------------------------------------------------------


pellipticalSlider =
function(type = c("persp", "contour"), B = 20)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively plots of probability
    
    # Description:
    #   Displays interactively plots of probability
    
    # Arguments:
    #   type - a character string specifying the plot type.
    #       Either a perspective plot which is the default or
    #       a contour plot with an underlying image plot will
    #       be created.
    #   B - the maximum slider menu value when the boundary
    #       value is infinite. By default this is set to 10.
    
    # FUNCTION:
    
    # Settings:
    type = match.arg(type)
    
    # Plot:
    if (type == "persp")
        .pellipticalPerspSlider(B = B)
    if (type == "contour")
        .pellipticalContourSlider(B = B)
        
    # Return Value:
    invisible()
}


# ------------------------------------------------------------------------------


.pnormCopula = 
function(u = 0.5, v = u, rho = 0.75, output = c("vector", "list") )
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes normal copula probability
    
    # Arguments:
    #   see function 'pellipticalCopula'
    
    # FUNCTION:
    
    # Type:
    output = match.arg(output)
    
    # Settings:
    type = "norm"
    if (is.list(u)) {
        v = u[[2]]
        u = u[[1]]
    }
    if (is.matrix(u)) {
        v = u[, 2]
        u = u[, 1]
    }
    
    # Copula Probability:
    C.uv = pnorm2d(qnorm(u), qnorm(v), rho = rho) 
    names(C.uv) = NULL
    
    # Simulates Max function:
    C.uv = (C.uv + abs(C.uv))/2       
    
    # On Boundary:
    C.uv[is.na(C.uv)] = 0      
    C.uv[which(u == 0)] = 0
    C.uv[which(u == 1)] = v[which(u == 1)]
    C.uv[which(v == 0)] = 0
    C.uv[which(v == 1)] = u[which(v == 1)]
    C.uv[which(u*v == 1)] = 1
    C.uv[which(u+v == 0)] = 0
    
    # Result:
    attr(C.uv, "control") <- c(rho = rho)
    
    # As List ?
    if (output == "list") {
        N = sqrt(length(u))
        x = u[1:N]
        y = matrix(v, ncol = N)[1, ]
        C.uv = list(x = x, y = y,  z = matrix(C.uv, ncol = N))
    }
    
    # Return Value:
    C.uv
}


# ------------------------------------------------------------------------------


.pcauchyCopula = 
function(u = 0.5, v = u, rho = 0.75, output = c("vector", "list") )
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes Student-t copula probability
    
    # Arguments:
    #   see function 'pellipticalCopula'
    
    # FUNCTION:
    
    # Cauchy Probability:
    C.uv = .ptCopula(u = u, v = v, rho = rho, nu = 1, output = output)
    attr(C.uv, "control") <- c(rho = rho)
    
    # Return Value:
    C.uv
}


# ------------------------------------------------------------------------------


.ptCopula = 
function(u = 0.5, v = u, rho = 0.75, nu = 4, output = c("vector", "list") )
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes Student-t copula probability
    
    # Arguments:
    #   see function 'pellipticalCopula'
    
    # FUNCTION:
    
    # Match Arguments:
    output = match.arg(output)
    
    # Settings:
    type = "t"
    if (is.list(u)) {
        v = u[[2]]
        u = u[[1]]
    }
    if (is.matrix(u)) {
        v = u[, 2]
        u = u[, 1]
    }
    
    # Copula Probability:
    C.uv = pt2d(qt(u, df = nu), qt(v, df = nu), rho = rho, nu = nu) 
    names(C.uv) = NULL
    
    # Simulates Max function:
    C.uv = (C.uv + abs(C.uv))/2       
    
    # On Boundary:
    C.uv[is.na(C.uv)] = 0      
    C.uv[which(u == 0)] = 0
    C.uv[which(u == 1)] = v[which(u == 1)]
    C.uv[which(v == 0)] = 0
    C.uv[which(v == 1)] = u[which(v == 1)]
    C.uv[which(u*v == 1)] = 1
    C.uv[which(u+v == 0)] = 0
    
    # Result:
    attr(C.uv, "control") <- c(rho = rho, nu = nu)
    
    # As List ?
    if (output == "list") {
        N = sqrt(length(u))
        x = u[1:N]
        y = matrix(v, ncol = N)[1, ]
        C.uv = list(x = x, y = y,  z = matrix(C.uv, ncol = N))
    }
    
    # Return Value:
    C.uv
}


# ------------------------------------------------------------------------------


.pellipticalCopulaGrid =
function(N, rho = 0.75, param = NULL, type = ellipticalList(), border = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes elliptical copula probability on a 2d grid
    
    # Arguments:
    #   see function pellipticalCopula()
    
    # FUNCTION:
    
    # Settings:
    U = (0:N)/N
    V = (1:(N-1))/N
    
    # Compute Density on Regular Grid:
    c.uv = .dellipticalCopulaGrid(N, rho, param, type, border = TRUE)
    c.uv$z[is.na(c.uv$z)] = 0      
    
    # Integrate to get Probability:
    if (TRUE) {
        C.uv = 0*c.uv$z
        for (i in 1:(N+1)) {
            for (j in 1:i) {    
                C.uv[i,j] = C.uv[j,i] = sum(c.uv$z[1:i, 1:j])
            }
        }
        C.uv = C.uv/N^2
    }
        
    if (FALSE) {
        # This is much slower !
        IJ = grid2d(1:(N+1))
        X = cbind(IJ$x, IJ$y)
        fun = function(X, C) sum(C[1:X[1], 1:X[2]])
        C.uv = apply(X, MARGIN=1, FUN = fun, C = c.uv$z)
        C.uv = matrix(C.uv, byrow = TRUE, ncol = N+1) / N^2
    } 
    
    # Probability - Take care about the Boundary on the Unit Square:
    C.uv[1, ] = C.uv[, 1] = 0
    C.uv[N+1, ] = C.uv[, N+1] = c.uv$x   
    names(C.uv) = NULL
    attr(C.uv, "control") <- c(rho = rho)
    C.uv = list(x = U, y = U, z = matrix(C.uv, ncol = length(U)))
    if (!border) {
        C.uv$z = C.uv$z[-1, ]
        C.uv$z = C.uv$z[-N, ]
        C.uv$z = C.uv$z[, -1]
        C.uv$z = C.uv$z[, -N]
        C.uv$x = C.uv$y = V
    }

    # Return Value:
    C.uv 
}


# ------------------------------------------------------------------------------


.pellipticalCopulaDiag =
function(N, rho = 0.75, param = NULL, type = ellipticalList(), border = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes elliptical diagonal cross section copula probability 
    
    # Arguments:
    #   see function pellipticalCopula()
    
    # FUNCTION:
    
    # Settings:
    U = (0:N)/N
    V = (1:(N-1))/N
    
    # Compute Density on Regular Grid:
    c.uv = .dellipticalCopulaGrid(N, rho, param, type[1], border = TRUE)
    c.uv$z[is.na(c.uv$z)] = 0      
    
    # Integrate to get Probability:
    C.uu = 0*U
    for (i in 1:(N+1)) {
        C.uu[i] = sum(c.uv$z[1:i, 1:i])
    }
    C.uu = C.uu/N^2
    names(C.uu) = NULL
    attr(C.uu, "control") <- c(rho = rho)
    
    if (border) {
        C.uu = list(x = U, y = C.uu)
    } else {
        C.uu = list(x = V, y = C.uu[c(-1,-(N+1))])
    }

    # Return Value:
    C.uu
}


# ------------------------------------------------------------------------------


.pellipticalPerspSlider =
function(B = 20)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively perspective plots of probability
    
    # Arguments:
    
    # FUNCTION:
    
    # Graphic Frame:
    par(mfrow = c(1, 1))
    
    # Internal Function:
    refresh.code = function(...)
    {
        # Startup Counter:
        .counter <<- .counter + 1
        if (.counter < 7) return ()
        
        # Sliders:
        Copula = .sliderMenu(no = 1)
        N = .sliderMenu(no = 2)
        rho = .sliderMenu(no = 3)
        nu = .sliderMenu(no = 4)
        s = .sliderMenu(no = 5)
        theta = .sliderMenu(no = 6)
        phi = .sliderMenu(no = 7)
        r = 1
        
        # Title:
        Names = 
            c("- Normal", "- Student t", "- Logistic", "- Exponential Power")   
        if (nu == 1) Names[2] = "- Student-t [Cauchy]"  
        if (s == 0.5) Names[4] = "- Exponential Power [Laplace]" 
        if (s == 1) Names[4] = "- Exponential Power [Kotz|Normal]"
        Title = paste("Elliptical Copula No:", as.character(Copula), 
            Names[Copula], "\nrho = ", as.character(rho)) 
        if (Copula == 2) Title = paste(Title, "nu =", as.character(nu))
        if (Copula == 4) Title = paste(Title, "s =", as.character(s))
        
        # Plot: 
        Type = c("norm", "t", "logistic", "epower")
        param = NULL
        if (Copula == 2) param = nu
        if (Copula == 4) param = c(r, s)
        P = .pellipticalCopulaGrid(N = N, rho = rho, param = param, 
            type = Type[Copula], border = TRUE)
        persp(P, theta = theta, phi = phi, col = "steelblue", shade = 0.5,
            ticktype = "detailed", cex = 0.5, xlab = "u", ylab = "v",
            zlab = "C(u, v)", xlim = c(0, 1), ylim = c(0, 1), zlim = c(0, 1) )
        title(main = Title)
        Tau = as.character(round(2*asin(rho)/pi, 2))
        mTitle = paste("Tau", Tau)
        mtext(mTitle, side = 4, col = "grey", cex = 0.7)
        mTitle = paste("1: Normal | 2: Student-t [Cauchy] | 3: Logistic |", 
            "4: Exponential Power [Laplace|Kotz]")
        mtext(mTitle, side = 1, line = 3, col = "grey", cex = 0.7)
                           
        # Reset Frame:
        par(mfrow = c(1, 1))
    }
  
    # Open Slider Menu:
    .counter <<- 0
    plot.names = c("Plot - theta", "... phi")
    .sliderMenu(refresh.code,
        names = c("Copula",  "N", "rho", "2: nu", "4: s", plot.names),
        minima      = c( 1,   10, -0.95,       1,    0.1,  -180,    0),
        maxima      = c( 4,  100,  0.95,       B,      5,   180,  360),
        resolutions = c( 1,   10,  0.05,     0.1,    0.1,     1,    1),
        starts      = c( 1,   20,  0.50,       4,      1,   -40,   30)) 
}


# ------------------------------------------------------------------------------


.pellipticalContourSlider =
function(B = 20)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively perspective plots of probability
    
    # Arguments:
    
    # FUNCTION:
    
    # Graphic Frame:
    par(mfrow = c(1, 1))
    
    # Internal Function:
    refresh.code = function(...)
    {
        # Startup Counter:
        .counter <<- .counter + 1
        if (.counter < 7) return ()
        
        # Sliders:
        Copula = .sliderMenu(no = 1)
        N = .sliderMenu(no = 2)
        rho = .sliderMenu(no = 3)
        nu = .sliderMenu(no = 4)
        s = .sliderMenu(no = 5)
        nlev = .sliderMenu(no = 6)
        ncol = .sliderMenu(no = 7)
        r = 1
        
        # Title:
        Names = 
            c("- Normal", "- Student t", "- Logistic", "- Exponential Power")   
        if (nu == 1) Names[2] = "- Student-t [Cauchy]"  
        if (s == 0.5) Names[4] = "- Exponential Power [Laplace]" 
        if (s == 1) Names[4] = "- Exponential Power [Kotz|Normal]"
        Title = paste("Elliptical Copula No:", as.character(Copula), 
            Names[Copula], "\nrho = ", as.character(rho)) 
        if (Copula == 2) Title = paste(Title, "nu =", as.character(nu))
        if (Copula == 4) Title = paste(Title, "s =", as.character(s))
        
        # Plot: 
        Type = c("norm", "t", "logistic", "epower")
        param = NULL
        if (Copula == 2) param = nu
        if (Copula == 4) param = c(r, s)
        P = .pellipticalCopulaGrid(N = N, rho = rho, param = param, 
            type = Type[Copula], border = FALSE)
        image(P, col = heat.colors(ncol), ylab = "v")
        mtext("u", side = 1, line = 2, cex = 0.7)
        contour(P, nlevels = nlev, add = TRUE)
        title(main = Title)
        Tau = as.character(round(2*asin(rho)/pi, 2))
        mTitle = paste("Tau", Tau)
        mtext(mTitle, side = 4, col = "grey", cex = 0.7)
        mTitle = paste("1: Normal | 2: Student-t [Cauchy] | 3: Logistic |", 
            "4: Exponential Power [Laplace|Kotz]")
        mtext(mTitle, side = 1, line = 3, col = "grey", cex = 0.7)
        
        # Reset Frame:
        par(mfrow = c(1, 1))
    }
  
    # Open Slider Menu:
    .counter <<- 0
    plot.names = c("Plot - levels", "... colors")
    .sliderMenu(refresh.code,
        names = c("Copula", "N", "rho", "2: nu", "4: s", plot.names),
        minima      = c( 1,  10, -0.95,       1,    0.1,     5,   12),
        maxima      = c( 4, 100,  0.95,       B,      5,   100,  256),
        resolutions = c( 1,  10,  0.05,     0.1,    0.1,     5,    4),
        starts      = c( 1,  20,  0.50,       4,      1,    10,   32)) 
}


################################################################################
# FUNCTION:                  ELLIPTICAL COPULAE DENSITY:
#  dellipticalCopula          Computes elliptical copula density 
#  dellipticalSlider          Generates interactive plots of density
#  .dnormCopula               Computes normal copula density
#  .dcauchyCopula             Computes Cauchy copula density
#  .dtCopula                  Computes Student-t copula density
#  .dellipticalCopulaGrid     Fast grid version for elliptical copula density
#  .dellipticalPerspSlider    Interactive perspective plots of density
#  .dellipticalContourSlider  Interactive contour plots of density


dellipticalCopula = 
function(u = 0.5, v = u, rho = 0.75, param = NULL, type = ellipticalList(), 
output = c("vector", "list"), border = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes extreme value copula density
    
    # Arguments:
    #   u, v - two numeric values or vectors of the same length at
    #       which the copula will be computed. If 'u' is a list then the
    #       the '$x' and '$y' elements will be used as 'u' and 'v'.
    #       If 'u' is a two column matrix then the first column will
    #       be used as 'u' and the the second as 'v'.
    #   rho - a numeric value setting the coorelation strength, ranging
    #       between minus one and one.
    #   param - additional distributional parameters.
    #   type - the type of the elliptical copula. Either "norm" or
    #       "t" denoting the normal or Student-t copula, respectively.
    #   output - a character string specifying how the output should
    #       be formatted. By default a vector of the same length as 
    #       'u' and 'v'. If specified as "list" then 'u' and 'v' are
    #       expected to span a two-dimensional grid as outputted by the 
    #       function 'grid2d' and the function returns a list with
    #       elements '$x', 'y', and 'z' which can be directly used 
    #       for example by 2D plotting functions. 
    
    # Value:
    #   returns a vector or list of probabilities depending on the
    #   value of the "output" variable.
    
    # Example:
    #   Diagonal Value: pnormCopula((0:10)/10)
    #   persp(pnormCopula(u = grid2d(), output = "list"))
    
    # FUNCTION:
    
    # Use Grid Version?
    if (is.numeric(u)) {
        if (length(u) == 1 & u[1] > 1) {
            ans = .dellipticalCopulaGrid(N = u, rho = rho, param = param, 
                type = type, border = border)
            return(ans)
        }
    }

    # Match Arguments:
    type = match.arg(type)
    output = match.arg(output)
    
    # Settings:
    if (is.list(u)) {
        v = u[[2]]
        u = u[[1]]
    }
    if (is.matrix(u)) {
        v = u[, 2]
        u = u[, 1]
    }
    if (length(u) == 1 & u[1] > 1) {
        return(.pellipticalCopulaGrid(N = u, rho, param, type, border = border))
    }
    
    # Parameters:
    if (type == "t") if (is.null(param)) param = c(nu = 4)
    if (type == "kotz") if (is.null(param)) param = c(r = 1)
    if (type == "epower") if (is.null(param)) param = c(r = 1, s = 1)
    
    # Density:  
    x = .qelliptical(u, param = param, type = type)
    y = .qelliptical(v, param = param, type = type)
    c.uv = delliptical2d(x, y, rho = rho, param = param, type = type) / (
        .delliptical(x, param = param, type = type) * 
        .delliptical(y, param = param, type = type) )
    if (rho == 0 & type == "norm") c.uv[!is.na(c.uv)] = 1
    names(c.uv) = NULL
    attr(c.uv, "control") <- c(rho = rho)
    if (output == "list") {
        N = sqrt(length(u))
        x = u[1:N]
        y = matrix(v, ncol = N)[1, ]
        c.uv = list(x = x, y = y, z = matrix(c.uv, ncol = N))
    }
    
    # Return Value:
    c.uv
}


# ------------------------------------------------------------------------------


dellipticalSlider =
function(type = c("persp", "contour"), B = 20)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively plots of density
    
    # Description:
    #   Displays interactively plots of density
    
    # Arguments:
    #   type - a character string specifying the plot type.
    #       Either a perspective plot which is the default or
    #       a contour plot with an underlying image plot will
    #       be created.
    #   B - the maximum slider menu value when the boundary
    #       value is infinite. By default this is set to 10.
    
    # FUNCTION:
    
    # Settings:
    type = match.arg(type)
    
    # Plot:
    if (type == "persp")
        .dellipticalPerspSlider(B = B)
    if (type == "contour")
        .dellipticalContourSlider(B = B)
        
    # Return Value:
    invisible()
}


# ------------------------------------------------------------------------------


.dnormCopula = 
function(u = 0.5, v = u, rho = 0.75, output = c("vector", "list") )
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes normal copula density
    
    # Arguments:
    #   see function 'dellipticalCopula'
    
    # FUNCTION:
    
    # Type:
    output = match.arg(output)
    
    # Settings:
    type = "norm"
    if (is.list(u)) {
        v = u[[2]]
        u = u[[1]]
    }
    if (is.matrix(u)) {
        v = u[, 2]
        u = u[, 1]
    }
    
    # Copula Density:
    x = qnorm(u)
    y = qnorm(v)
    c.uv = dnorm2d(x, y, rho)/(dnorm(x) * dnorm(y))
    names(c.uv) = NULL
    
    # Result:
    attr(c.uv, "control") <- c(rho = rho)
    
    # As List ?
    if (output == "list") {
        N = sqrt(length(u))
        x = u[1:N]
        y = matrix(v, ncol = N)[1, ]
        c.uv = list(x = x, y = y,  z = matrix(c.uv, ncol = N))
    }
    
    # Return Value:
    c.uv
}


# ------------------------------------------------------------------------------


.dtCopula = 
function(u = 0.5, v = u, rho = 0.75, nu = 4, output = c("vector", "list") )
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes Student-t copula density
    
    # Arguments:
    #   see function 'dellipticalCopula'
    
    # FUNCTION:
    
    # Match Arguments:
    output = match.arg(output)
    
    # Settings:
    type = "t"
    if (is.list(u)) {
        v = u[[2]]
        u = u[[1]]
    }
    if (is.matrix(u)) {
        v = u[, 2]
        u = u[, 1]
    }
    
    # Copula Probability:
    x = qt(u, df = nu)
    y = qt(v, df = nu)
    c.uv = dt2d(x, y, rho, nu)/(dt(x, nu) * dt(y, nu))
    names(c.uv) = NULL
    
    # Result:
    attr(c.uv, "control") <- c(rho = rho, nu = nu)
    
    # As List ?
    if (output == "list") {
        N = sqrt(length(u))
        x = u[1:N]
        y = matrix(v, ncol = N)[1, ]
        c.uv = list(x = x, y = y,  z = matrix(c.uv, ncol = N))
    }
    
    # Return Value:
    c.uv
}


# ------------------------------------------------------------------------------


.dcauchyCopula = 
function(u = 0.5, v = u, rho = 0.75, nu = 4, output = c("vector", "list") )
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes Student-t copula density
    
    # Arguments:
    #   see function 'dellipticalCopula'
    
    # FUNCTION:
    
    # Cauchy Density:
    c.uv = .dtCopula(u = u, v = v, rho = rho, nu = 1, output = output)
    attr(c.uv, "control") <- c(rho = rho)
    
    # Return Value:
    c.uv
}


# ------------------------------------------------------------------------------


.dellipticalCopulaGrid =
function(N, rho = 0.75, param = NULL, type = ellipticalList(), border = TRUE)
{   # A function implemented by Diethelm Wuertz

    # Description:
    #   Computes extreme value copula density
    
    # Arguments:
    #   N - the number of grid points is (N+1)*(N+1)
    #   rho - a numeric value setting the coorelation strength, ranging
    #       between minus one and one.
    #   param - additional distributional parameters.
    #   type - the type of the elliptical copula. Either "norm" or
    #       "t" denoting the normal or Student-t copula, respectively.
    
    # Value:
    #   returns a vector or list of probabilities depending on the
    #   value of the "output" variable.
    
    # Note:
    #   Made for the Sliders.
    
    # FUNCTION:
    
    # Settings:
    type = type[1]
    U = (0:N)/N
    V = (1:(N-1))/N

    # Reduce to Grid - speeds up the computation:
    M = N%/%2 + 1
    X = .qelliptical(U[1:M], param = param, type = type)
    if (N%%2 == 0) {
        X = c(X, rev(-X)[-1]) 
    } else {
        X = c(X, rev(-X)) 
    }
    NX = length(X)
    x = rep(X, times = NX)
    y = rep(X, each = NX)
    D = .delliptical(X, param = param, type = type)
    DX = rep(D, times = NX)
    DY = rep(D, each = NX)

    # Density:  
    c.uv = delliptical2d(x, y, rho = rho, param = param, type = type) / (DX*DY)
    if (rho == 0 & type == "norm") c.uv[!is.na(c.uv)] = 1
    c.uv[is.na(c.uv)] = 0
    names(c.uv) = NULL
    attr(c.uv, "control") <- c(rho = rho)
    c.uv = list(x = U, y = U, z = matrix(c.uv, ncol = N+1))
    if (!border) {
        c.uv$z = c.uv$z[-1, ]
        c.uv$z = c.uv$z[-N, ]
        c.uv$z = c.uv$z[, -1]
        c.uv$z = c.uv$z[, -N]
        c.uv = list(x = V, y = V, z = matrix(c.uv$z, ncol = N-1))
    }
    
    # Return Value:
    c.uv
}


# ------------------------------------------------------------------------------


.dellipticalPerspSlider =
function(B = 20)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively perspective plots of density
    
    # FUNCTION:
    
    # Graphic Frame:
    par(mfrow = c(1, 1))
    
    # Internal Function:
    refresh.code = function(...)
    {
        # Startup Counter:
        .counter <<- .counter + 1
        if (.counter < 7) return ()
        
        # Sliders:
        Copula = .sliderMenu(no = 1)
        N = .sliderMenu(no = 2)
        rho = .sliderMenu(no = 3)
        nu = .sliderMenu(no = 4)
        s = .sliderMenu(no = 5)
        theta = .sliderMenu(no = 6)
        phi = .sliderMenu(no = 7)
        r = 1
        
        # Title:
        Names = 
            c("- Normal", "- Student t", "- Logistic", "- Exponential Power")   
        if (nu == 1) Names[2] = "- Student-t [Cauchy]"  
        if (s == 0.5) Names[4] = "- Exponential Power [Laplace]" 
        if (s == 1) Names[4] = "- Exponential Power [Kotz|Normal]"
        Title = paste("Elliptical Copula Density No:", as.character(Copula), 
            Names[Copula], "\nrho = ", as.character(rho)) 
        if (Copula == 2) Title = paste(Title, "nu =", as.character(nu))
        if (Copula == 4) Title = paste(Title, "s =", as.character(s))
        
        # Plot: 
        uv = grid2d(x = (1:(N-1))/N)
        Type = c("norm", "t", "logistic", "epower")
        param = NULL
        if (Copula == 2) param = nu
        if (Copula == 4) param = c(r, s)
        D = .dellipticalCopulaGrid(N, rho = rho, param = param, 
            type = Type[Copula], border = FALSE)
        Integrated = as.character(round(mean(D$z),2))
        Var = var(as.vector(D$z), na.rm = TRUE)
        if (Var < 1.0e-6) {
            # A flat perspective plot fails, if zlim is not specified!
            Mean = round(1.5*mean(as.vector(D$z), na.rm = TRUE), 2)
            persp(D, theta = theta, phi = phi, col = "steelblue", shade = 0.5,
                ticktype = "detailed", cex = 0.5, xlab = "u", ylab = "v",
                zlim = c(0, Mean), zlab = "C(u,v)" )
        } else {
            persp(D, theta = theta, phi = phi, col = "steelblue", shade = 0.5,
                ticktype = "detailed", cex = 0.5, xlab = "u", ylab = "v",
                zlab = "C(u,v)" )
        }
        title(main = Title)
        Tau = as.character(round(2*asin(rho)/pi, 2))
        mTitle = paste("Mean: ", Integrated, " |  Tau", Tau)
        mtext(mTitle, side = 4, col = "grey", cex = 0.7)
        mTitle = paste("1: Normal | 2: Student-t [Cauchy] | 3: Logistic |", 
            "4: Exponential Power [Laplace|Kotz]")
        mtext(mTitle, side = 1, col = "grey", cex = 0.7)
                           
        # Reset Frame:
        par(mfrow = c(1, 1))
    }
  
    # Open Slider Menu:
    .counter <<- 0
    plot.names = c("Plot - theta", "... phi")
    .sliderMenu(refresh.code,
        names = c("Copula", "N", "rho", "3: nu", "4: s", plot.names),
        minima      = c( 1,  10, -0.95,       1,   0.1,  -180,    0),
        maxima      = c( 4, 100,  0.95,       B,     5,   180,  360),
        resolutions = c( 1,  10,  0.05,     0.1,   0.1,     1,    1),
        starts      = c( 1,  20,  0.50,       4,     1,   -40,   30)) 
}


# ------------------------------------------------------------------------------


.dellipticalContourSlider =
function(B = 20)
{   # A function implemented by Diethelm Wuertz
        
    # Description:
    #   Displays interactively perspective plots of density
    
    #FUNCTION:
    
    # Graphic Frame:
    par(mfrow = c(1, 1))
    
    # Internal Function:
    refresh.code = function(...)
    {
        # Startup Counter:
        .counter <<- .counter + 1
        if (.counter < 7) return ()
        
        # Sliders:
        Copula = .sliderMenu(no = 1)
        N = .sliderMenu(no = 2)
        rho = .sliderMenu(no = 3)
        nu = .sliderMenu(no = 4)
        s = .sliderMenu(no = 5)
        nlev = .sliderMenu(no = 6)
        ncol = .sliderMenu(no = 7)
        if (rho == 0 & Copula == 1) return(invisible())
        r = 1
        
        # Title:
        Names = 
            c("- Normal", "- Student t", "- Logistic", "- Exponential Power")   
        if (nu == 1) Names[2] = "- Student-t [Cauchy]"  
        if (s == 0.5) Names[4] = "- Exponential Power [Laplace]" 
        if (s == 1) Names[4] = "- Exponential Power [Kotz|Normal]"
        Title = paste("Elliptical Copula Density No:", as.character(Copula), 
            Names[Copula], "\nrho = ", as.character(rho)) 
        if (Copula == 2) Title = paste(Title, "nu =", as.character(nu))
        if (Copula == 4) Title = paste(Title, "s =", as.character(s))
        
        # Plot: 
        uv = grid2d(x = (0:N)/N)
        Type = c("norm", "t", "logistic", "laplace", "kotz", "epower")
        param = NULL
        if (Copula == 2) param = nu
        if (Copula == 5) param = r
        if (Copula == 6) param = c(r, s)
        D = .dellipticalCopulaGrid(N, rho = rho, param = param, 
            type = Type[Copula], border = FALSE)
        Integrated = as.character(round(mean(D$z),2))
        image(D, col = heat.colors(ncol), ylab = "v",
            xlim = c(0,1), ylim = c(0,1) )
        mtext("u", side = 1, line = 2, cex = 0.7)
        contour(D, nlevels = nlev, add = TRUE)
        title(main = Title)
        Tau = as.character(round(2*asin(rho)/pi, 2))
        mTitle = paste("Mean: ", Integrated, " |  Tau", Tau)
        mtext(mTitle, side = 4, col = "grey", cex = 0.7)
        mTitle = paste("1: Normal | 2: Student-t [Cauchy] | 3: Logistic |", 
            "4: Exponential Power [Laplace|Kotz]")
        mtext(mTitle, side = 1, line =  3, col = "grey", cex = 0.7)
        
        # Reset Frame:
        par(mfrow = c(1, 1))
    }
  
    # Open Slider Menu:
    .counter <<- 0
    plot.names = c("Plot - levels", "... colors")
    .sliderMenu(refresh.code,
        names = c("Copula", "N", "rho", "2: nu", "4: s", plot.names),
        minima      = c( 1,  10, -0.95,       1,    0.1,     5,   12),
        maxima      = c( 4, 100,  0.95,       B,      5,   100,  256),
        resolutions = c( 1,  10,  0.05,     0.1,    0.1,     5,    4),
        starts      = c( 1,  20,  0.50,       4,      1,    10,   32)) 
}


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