The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Densities of Two-Level Nested Archimedean Copulas

Marius Hofert

2025-02-09

require(copula)
require(grid)
require(lattice)
source(system.file("Rsource", "dnac.R", package="copula"))
set.seed(271)

Note that nacLL() will use package partitions and polynom.

Examples (sampling and evaluating the log-likelihood)

Example 1: ((1,2), (3,4,5))-Gumbel

Consider the following setup.

n <- 250
family <- "Gumbel"
tau <- c(0.2, 0.4, 0.6)
th <- getAcop(family)@iTau(tau)
G3 <- onacopula(family, C(th[1], , list(C(th[2], 1:2), C(th[3], 3:5))))

Sample and compute the log-likelihood.

U <- rnacopula(n, G3)
nacLL(G3, u=U) # log-likelihood at correct parameters
## [1] 450.654

Example 2: (1, (2,3), 4, (5,6,7))-Gumbel

Consider the following setup.

n <- 250
family <- "Gumbel"
cop. <- getAcop(family)
tau <- c(0.2, 0.4, 0.6)
th <- cop.@iTau(tau)
cop <- onacopula(family, C(th[1], c(1,4),
                           list(C(th[2], 2:3), C(th[3], 5:7))))

Sample and compute the log-likelihood.

U <- rnacopula(n, cop)
nacLL(cop, u=U) # log-likelihood at correct parameters
## [1] 455.4679

Example 3: (1, (2,3))-Gumbel

Consider the following setup.

n <- 250
family <- "Gumbel"
tau <- c(0.25, 0.5)
th <- getAcop(family)@iTau(tau)
copTrue <- onacopula(family, C(th[1], 1, C(th[2], 2:3)))

Sample and compute the log-likelihood.

U <- rnacopula(n, copTrue)
nacLL(copTrue, u=U) # log-likelihood at correct parameters
## [1] 136.6138

Plots of the negative log-likelihood

We consider a (1, (2,..,\(d\)))-structure structure here (we choose \(d=3\) here but larger \(d\) are of course possible; plotting can be done as long as we consider two parameters only).

family <- "Gumbel" # choose "Clayton" or "Gumbel"
 compTr <- 1 # non-sectorial indices; *Tr stands for the true (nesting structure/model)
scompTr <- 2:3 # sectorial indices (for plotting, need 2:d)
stopifnot(compTr==1) # otherwise, sub (for plotting, see below) is wrong

We first define the negative log-likelihood.

##' Negative Log Likelihood for the two-parameter case
##' C_0({u_j}, C_1({u_k})) where j in 'comp';  k in 'scomp'
nLL2 <- function(th, u, family, comp, scomp)
{
    stopifnot(length(th) == 2)
    if(th[1] > th[2]) # sufficient nesting condition not fulfilled
    return(Inf) # for minimization
    cop <- onacopulaL(family, list(th[1], comp, list(list(th[2], scomp))))
    -nacLL(cop, u=u)
}

Determine the values of the negative log-likelihood on a grid

n <- 100
cop. <- getAcop(family)
tau <- c(0.25, 0.5)
(thTr <- cop.@iTau(tau))
## [1] 1.333333 2.000000
cop <- onacopula(family, C(thTr[1], compTr, C(thTr[2], scompTr))) # copula
U <- rnacopula(n, cop) # sample
h <- 0.2 # delta{tau} for defining a range of theta's
(th0 <- cop.@iTau(c(tau[1]-h, tau[1]+h)))
## [1] 1.052632 1.818182
(th1 <- cop.@iTau(c(tau[2]-h, tau[2]+h)))
## [1] 1.428571 3.333333
m <- 20 # number of grid points
grid <- expand.grid(th0= seq(th0[1], th0[2], length.out=m),
                    th1= seq(th1[1], th1[2], length.out=m))
val.grid <- apply(grid, 1, nLL2,
          u=U, family=family, comp=compTr, scomp=scompTr)

Plotting

First we determine some plot supplements including the optimum of the negative log-likelihood on the grid.

true.val <- c(th0=thTr[1], th1=thTr[2],
              nLL=nLL2(thTr, u=U, family=family, comp=compTr, scomp=scompTr)) # true value
ind <- which.min(val.grid)
opt.grd <- c(grid[ind,], nLL=val.grid[ind]) # optimum on the grid
pts <- rbind(true.val, opt.grd) # points to add to wireframe plot
title <- paste("-log-likelihood of a nested", family, "copula") # title
mysec <- { if(length(scompTr)==2) bquote(italic(u[3]))
           else substitute(list(...,italic(u[j])), list(j=max(scompTr))) }
sub <- substitute(italic(C(bolditalic(u)))==italic(C[0](u[1],C[1](u[2],MSEC))) ~~~~~~
                  italic(n)==N ~~~~~~ tau(theta[0])==TAU0 ~~~~~~ tau(theta[1])==TAU1,
                  list(MSEC=mysec, N=n, TAU0=tau[1], TAU1=tau[2]))
sub <- as.expression(sub) # lattice "bug" (only needed by lattice)
xlab <- expression(italic(theta[0]))
ylab <- expression(italic(theta[1]))
zlab <- list(as.expression(-log~L *
    group("(",italic(theta[0])*"," ~ italic(theta[1])*";"~bolditalic(u),")")),
             rot = 90)
sTit <- list(c(expression(group("(",list(theta[0],theta[1]),")")^T),
           expression(group("(",list(hat(theta)["0,n"],hat(theta)["1,n"]),")")^T)))

Now consider a wireframe and a level plot.

wireframe(val.grid~grid[,1]*grid[,2], aspect=1, zoom=1.02, xlim=th0, ylim=th1,
          zlim= range(val.grid, as.numeric(pts[,3]), finite=TRUE),
          xlab=xlab, ylab=ylab, zlab=zlab, main=title, sub=sub, pts=pts,
          par.settings=list(standard.theme(color=FALSE),
          layout.heights=list(sub=2.4), background=list(col="#ffffff00"),
          axis.line=list(col="transparent"), clip=list(panel="off")),
          alpha.regions=0.5, scales=list(col=1, arrows=FALSE),
          ## add wire/points
          panel.3d.wireframe = function(x, y, z, xlim, ylim, zlim, xlim.scaled,
          ylim.scaled, zlim.scaled, pts, ...) {
              panel.3dwire(x=x, y=y, z=z, xlim=xlim, ylim=ylim, zlim=zlim,
                           xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
                           zlim.scaled=zlim.scaled, ...)
              panel.3dscatter(x=as.numeric(pts[,1]), y=as.numeric(pts[,2]),
                              z=as.numeric(pts[,3]), xlim=xlim, ylim=ylim, zlim=zlim,
                              xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
                              zlim.scaled=zlim.scaled, type="p", col=1,
                              pch=c(3,4), lex=2, cex=1.4, .scale=TRUE, ...)
          },
          key =
              list(x=-0.01, y=1, points=list(pch=c(3,4), col=1, lwd=2, cex=1.4),
                   text = sTit, padding.text=3, cex=1, align=TRUE, transparent=TRUE))

levelplot(val.grid~grid[,1]*grid[,2], aspect=1, xlab=xlab, ylab=ylab,
          par.settings=list(layout.heights=list(main=3, sub=2),
          regions=list(col=gray(140:400/400))),
          xlim= extendrange(grid[,1], f = 0.04),
          ylim= extendrange(grid[,2], f = 0.04),
          main=title, sub=sub, pts=pts,
          scales=list(alternating=c(1,1), tck=c(1,0)), contour=TRUE,
          panel=function(x, y, z, pts, ...){
              panel.levelplot(x=x, y=y, z=z, ...)
              grid.points(x=pts[1,1], y=pts[1,2], pch=3,
                          gp=gpar(lwd=2, col="black")) # + true value
              grid.points(x=pts[2,1], y=pts[2,2], pch=4,
                          gp=gpar(lwd=2, col="black")) # x optimum
          },
          key =
              list(x=0.18, y=1.09, points=list(pch=c(3,4), col=1, lwd=2, cex=1.4),
                   columns = 2, text = sTit, align=TRUE, transparent=TRUE))

Computing the MLE (via optimization)

For illustration purposes, we start not too closely to the optimum.

ropt <- optim(c(1, 3), nLL2,
              u=U, family=family, comp=compTr, scomp=scompTr)

Now compare the different results (true values, optimum on the grid, optimum via optim()).

rbind(pts, optim=c(ropt$par, ropt$value))
##          th0      th1      nLL      
## true.val 1.333333 2        -82.54895
## opt.grd  1.374969 1.929825 -82.81026
## optim    1.361803 1.915033 -82.82763

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.