Generalized Inverse Gaussian Archimedean Copulas

Marius Hofert

2019-04-21

source(system.file("Rsource", "GIG.R", package="copula"))## ../inst/Rsource/GIG.R
require(copula)
require(bbmle)
require(lattice)
require(grid)
doPDF <- FALSE ## set 'do.profile' below -- *visibly*

1 Auxiliary functions

##' Initial interval for GIG
##' @title Initial interval for GIG
##' @param U (n x d)-matrix of simulated data
##' @param h non-negative auxiliary parameter for computing initial intervals
##' @param method "etau" via sample version of Kendall's tau
##'               "dmle.G" via DMLE of Gumbel
##' @return (2 x 2)-matrix containing the initial interval [1st row: lower,
##'         2nd row: upper; 2 parameters => 2 cols]
##' @author Marius Hofert
ii.GIG <- function(U, h, method=c("etau","dmle.G")){
    stopifnot(h >= 0, length(h) >= 2)
    I <- matrix(, nrow=2, ncol=2, dimnames=list(c("lower", "upper"), c("nu", "theta")))
    ## estimate Kendall's tau
    method <- match.arg(method)
    tau.hat <- switch(method,
                      "etau" = { # uses sample version of tau, more accurate but slower
                          tau.hat.mat <- cor(U, method="kendall")
                          mean(tau.hat.mat[upper.tri(tau.hat.mat)])
                      },
                      "dmle.G" = { # uses DMLE for Gumbel to estimate tau
                          Z <- apply(U, 1, max)
                          theta.hat.G <- log(ncol(U))/(log(length(Z))-log(sum(-log(Z))))
                          copGumbel@tau(theta.hat.G)
                      },
                      stop("wrong method:", method))
    ## compute largest value of theta (for upper left endpoint of the inital interval)
    stopifnot(tau.hat > 0)
    nu.min <- 0
    I[1,1] <- nu.min # smallest value for nu
    th.max <- iTau.GIG(max(tau.hat-h[1],0.005), theta=c(nu.min, NA))
    I[2,2] <- th.max[2] # largest value for theta
    ## compute smallest theta (for lower left endpoint of the inital interval)
    th.min <- iTau.GIG(min(tau.hat+h[2],0.995), theta=c(nu.min, NA)) # largest attainable tau with 1e-30 is one.m.eps=0.9602
    I[1,2] <- th.min[2]
    ## compute largest nu (for lower right endpoint of the inital interval)
    nu.max <- iTau.GIG(max(tau.hat-h[1],0.005), theta=c(NA, th.min[2]))
    I[2,1] <- nu.max[1]
    ## result
    I
}
##' -log-likelihood
##' @title -log-likelihood
##' @param nu parameter of the generator/copula
##' @param theta parameter of the generator/copula
##' @param u (n x d)-matrix of simulated data
##' @return -sum(log(density))
##' @author Marius Hofert
nlogl.GIG <- function(nu, theta, u){
    if(!is.matrix(u)) u <- rbind(u)
    if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2
    -sum(dacopula.GIG(u, theta=c(nu, theta), n.MC=0, log=TRUE))
}
nlogl.GIG. <- function(theta, u) nlogl.GIG(theta[1], theta=theta[2], u=u) # vectorized version

2 Setup

Note: The GIG family is two-parametric.

Plot Kendall’s tau as a function in \(\theta\) for different \(\nu\)

Parameter specification

Let’s specify some parameters.

3 Sampling and estimation

Sampling

Estimation

##    user  system elapsed 
##  67.949   0.525  68.938
##    user  system elapsed 
##  49.728   0.027  50.063
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nLL, start = list(nu = mean(I[, 1]), theta = mean(I[, 
##     2])), method = "L-BFGS-B", lower = c(nu = I[1, 1], theta = I[1, 
##     2]), upper = c(nu = I[2, 1], theta = I[2, 2]))
## 
## Coefficients:
##         Estimate Std. Error
## nu    0.19074320 0.06197105
## theta 0.09566594 0.01386170
## 
## -2 log L: -994.5633
## List of 6
##  $ par        : Named num [1:2] 0.1907 0.0957
##   ..- attr(*, "names")= chr [1:2] "nu" "theta"
##  $ value      : num -497
##  $ counts     : Named int [1:2] 30 30
##   ..- attr(*, "names")= chr [1:2] "function" "gradient"
##  $ convergence: int 0
##  $ message    : chr "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
##  $ hessian    : num [1:2, 1:2] 264 139 139 5277
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "nu" "theta"
##   .. ..$ : chr [1:2] "nu" "theta"
##    user  system elapsed 
##  59.434   0.013  59.825
## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = nlogl.GIG, start = list(nu = mean(I[, 1]), theta = mean(I[, 
##     2])), method = "L-BFGS-B", data = list(u = U), lower = c(nu = I[1, 
##     1], theta = I[1, 2]), upper = c(nu = I[2, 1], theta = I[2, 
##     2]))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## nu    0.190743   0.060871  3.1336  0.001727 ** 
## theta 0.095666   0.013684  6.9910 2.729e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: -994.5633
## List of 8
##  $ par        : Named num [1:2] 0.1907 0.0957
##   ..- attr(*, "names")= chr [1:2] "nu" "theta"
##  $ value      : num -497
##  $ counts     : Named int [1:2] 30 30
##   ..- attr(*, "names")= chr [1:2] "function" "gradient"
##  $ convergence: int 0
##  $ message    : chr "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
##  $ hessian    : num [1:2, 1:2] 272 112 112 5387
##  $ maxgrad    : num 9.36
##  $ eratio     : num 0.0501

4 Plots

Profile likelihood plots

-log-likelihood plot

##    user  system elapsed 
## 190.988   0.049 192.290

Wireframe

true.theta <- theta
true.val <- c(true.theta, nlogl.GIG.(true.theta, u=U)) # theoretical optimum
opt <- ml@coef # optimizer-optimum
opt.val <- c(opt, nlogl.GIG.(opt, u=U)) # optimizer-optimum and its value
pts <- rbind(true.val, opt.val) # points to add to wireframe plot
title <- "-log-likelihood of an Archimedean GIG copula" # title
sub <- substitute(italic(n) == N ~~~  italic(d)== D ~~~
                  tau == TAU ~~~ "#{eval}:" ~ NIT,
                  list(N=n, D=d, TAU= tau, NIT= ml@details$counts[[1]]))
sub <- as.expression(sub) # lattice bug
wireframe(val.grid~grid[,1]*grid[,2], screen=list(z=70, x=-55), zoom=0.95,
          xlab=expression(italic(theta)), ylab=expression(italic(beta)),
          zlab = list(as.expression(-log~L * group("(",list(theta,beta),")")), rot=90),
          main=title, sub=sub, pts=pts, scales=list(col=1, arrows=FALSE),
          par.settings=list(axis.line=list(col="transparent"),
          clip=list(panel="off")), zlim=c(min(val.grid, pts[,3]),
                                   max(val.grid, pts[,3])), aspect=1,
          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, alpha.regions=0.8, ...)
              panel.3dscatter(x=pts[,1], y=pts[,2], z=pts[,3],
                              xlim=xlim, ylim=ylim, zlim=zlim,
                              xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
                              zlim.scaled=zlim.scaled, type="p", col=c("red","blue"),
                              pch=c(3,4), lex=2, cex=1.4, .scale=TRUE, ...)
          }, key=list(x=0.64, y=1.01, points=list(pch=c(3,4), col=c("red","blue"),
                                      lwd=2, cex=1.4),
             text=list(c("True value", "Optimum of optimizer")), padding.text=3,
             cex=1, align=TRUE, transparent=TRUE))

Levelplot