Densities of Two-Level Nested Archimedean Copulas

Marius Hofert

2018-12-20

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.

Sample and compute the log-likelihood.

## [1] 450.654

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

Consider the following setup.

Sample and compute the log-likelihood.

## [1] 455.4679

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

Consider the following setup.

Sample and compute the log-likelihood.

## [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

## [1] 1.333333 2.000000
## [1] 1.052632 1.818182
## [1] 1.428571 3.333333

Plotting

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

Now consider a wireframe and a level plot.

Computing the MLE (via optimization)

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

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

##          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