copula (version 0.999-7)

emle: Maximum Likelihood Estimators for (Nested) Archimedean Copulas

Description

Compute (simulated) maximum likelihood estimators for (nested) Archimedean copulas.

Usage

emle(u, cop, n.MC=0, optimizer="optimize", method,
     interval=initOpt(cop@copula@name),
     start=list(theta=initOpt(cop@copula@name, interval=FALSE, u=u)),
     ...)
.emle(u, cop, n.MC=0,
      interval=initOpt(cop@copula@name), ...)

Arguments

u

$n\times d$-matrix of (pseudo-)observations (each value in $[0,1]$) from the copula, with $n$ the sample size and $d$ the dimension.

cop

outer_nacopula to be estimated (currently only non-nested, that is, latex{Archi-medean}{Archimedean} copulas are admitted).

n.MC

integer, if positive, simulated maximum likelihood estimation (SMLE) is used with sample size equal to n.MC; otherwise (n.MC=0), MLE. In SMLE, the $d$th gener

optimizer

a string or NULL, indicating the optimizer to be used, where NULL means to use optim via the standard Rfunction mle() from

method

only when optimizer is NULL or "optim", the method to be used for optim.

interval

bivariate vector denoting the interval where optimization takes place. The default is computed as described in Hofert et al. (2011a).

start

list of initial values, passed through.

...

additional parameters passed to optimize.

Value

• [object Object],[object Object]

Details

Exact formulas for the generator derivatives were derived in Hofert et al. (2011b). Based on these formulas one can compute the (log-)densities of the Archimedean copulas. Note that for some densities, the formulas are numerically highly non-trivial to compute and considerable efforts were put in to make the computations numerically feasible even in large dimensions (see the source code of the Gumbel copula, for example). Both MLE and SMLE showed good performance in the simulation study conducted by Hofert et al. (2011a) including the challenging 100-dimensional case. Alternative estimators (see also enacopula) often used because of their numerical feasibility, might break down in much smaller dimensions.

Note: SMLE for Clayton currently faces serious numerical issues and is due to further research. This is only interesting from a theoretical point of view, since the exact derivatives are known and numerically non-critical to evaluate.

References

Hofert, M., Mächler{Maechler}, M., and McNeil, A. J. (2011a). Estimators for Archimedean copulas in high dimensions: A comparison; submitted.

Hofert, M., Mächler{Maechler}, M., and McNeil, A. J. (2012). Likelihood inference for Archimedean copulas in high dimensions under known margins. Journal of Multivariate Analysis 110, 133--150.

See Also

mle2 from package bbmle and mle from stats4 on which mle2 is modeled. enacopula (wrapper for different estimators). demo(opC-demo) and demo(GIG-demo) for examples of two-parameter families.

Examples

tau <- 0.25
(theta <- copGumbel@iTau(tau)) # 4/3
d <-  20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))

set.seed(1)
n <- 200
U <- rnacopula(n,cop)

## Estimation
system.time(efm <- emle(U, cop))
summary(efm)

## Profile likelihood plot
if(getRversion() >= "2.15.0") { # otherwise, need  'stats4::' everywhere
pfm <- profile(efm)
ci  <- confint(pfm, level=0.95)
print(ci)
stopifnot(ci[1] <= theta, theta <= ci[2])
plot(pfm)               # |z| against theta, |z| = sqrt(deviance)
plot(pfm, absVal=FALSE, #  z  against theta
      show.points=TRUE) # showing how it's interpolated
## and show the true theta:
abline(v=theta, col="lightgray", lwd=2, lty=2)
axis(1, pos = 0, at=theta, label=expression(theta[0]))

## Plot of the log-likelihood, MLE  and  conf.int.:
logL <- function(x) -efm@minuslogl(x)
       # == -sum(copGumbel@dacopula(U, theta=x, log=TRUE))
logL. <- Vectorize(logL)
I <- c(cop@copula@iTau(0.1), cop@copula@iTau(0.4))
curve(logL., from=I[1], to=I[2], xlab=quote(theta),
      ylab="log-likelihood",
      main="log-likelihood for Gumbel")
abline(v = c(theta, efm@coef), col="magenta", lwd=2, lty=2)
axis(1, at=c(theta, efm@coef), padj = c(-0.5, -0.8), hadj = -0.2,
     col.axis="magenta", label= expression(theta[0], hat(theta)[n]))
abline(v=ci, col="gray30", lwd=2, lty=3)
text(ci[2], extendrange(par("usr")[3:4], f= -.04)[1],
     "95% conf. int.", col="gray30", adj = -0.1)
}#-- only in R >= 2.15.0