# Densities of Two-Level Nested Archimedean Copulas

#### 2023-01-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.

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),
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