estim.misc {copula} | R Documentation |
Various Estimators for (Nested) Archimedean Copulas
Description
Various Estimators for (Nested) Archimedean Copulas, namely,
- ebeta
Method-of-moments-like estimator based on (a multivariate version of) Blomqvist'sbeta.
- edmle
Maximum likelihood estimator based on the diagonal of a (nested) Archimedean copula.
- etau
Method-of-moments-like estimators based on (bivariate) Kendall's tau.
Usage
ebeta(u, cop, interval = initOpt(cop@copula@name), ...)
edmle(u, cop, interval = initOpt(cop@copula@name), warn=TRUE, ...)
etau(u, cop, method = c("tau.mean", "theta.mean"), warn=TRUE, ...)
Arguments
u |
|
cop |
|
interval |
bivariate vector denoting the interval where optimization takes place. The default is computed as described in Hofert et al. (2013). |
method |
a character string specifying the method (only
for
|
warn |
logical indicating if warnings are printed:
|
... |
additional arguments passed to
|
Details
For ebeta
, the parameter is estimated with a
method-of-moments-like procedure such that the population version of
the multivariate Blomqvist's beta matches its sample version.
Note that the copula diagonal is a distribution function and the
maximum of all components of a random vector following the copula is
distributed according to this distribution function. For
edmle
, the parameter is estimated via maximum-likelihood
estimation based on the diagonal.
For etau
, corKendall(u, ...)
is used and if there
are no NA
s in u
, by default (if no additional
arguments are provided), corKendall()
calls the O(n log(n))
fast cor.fk()
from package pcaPP
instead of the O(n^2)
cor(*, method="kendall")
.
Conversely, when u
has NA
s, by default,
corKendall(u, ...)
will use
cor(u, method="kendall", use = "pairwise")
such that
etau(u, *)
will work.
Furthermore, method="tau.mean"
means that the average
of sample versions of Kendall's tau are computed first and then the
parameter is determined such that the population version of Kendall's
tau matches this average (if possible); the method="theta.mean"
stands for first computing all pairwise Kendall's tau estimators and
then returning the mean of these estimators.
For more details, see Hofert et al. (2013).
Note that these estimators should be used with care; see the
performance results in Hofert et al. (2013). In particular,
etau
should be used with the (default) method "tau.mean"
since "theta.mean"
is both slower and more prone to errors.
Value
ebeta
the return value of
safeUroot
(that is, typically almost the same as the value ofuniroot
) giving the Blomqvist beta estimator.edmle
list
as returned byoptimize
, including the diagonal maximum likelihood estimator.etau
method-of-moments-like estimator based on Kendall's tau for the chosen method.
References
Hofert, M., Mächler, M., and McNeil, A. J. (2013). Archimedean Copulas in High Dimensions: Estimators and Numerical Challenges Motivated by Financial Applications. Journal de la Société Française de Statistique 154(1), 25–63.
See Also
corKendall()
.
The more sophisticated estimators emle
(Maximum Likelihood) and
emde
(Minimum Distance). enacopula
(wrapper for different estimators).
Examples
tau <- 0.25
(theta <- copGumbel@iTau(tau)) # 4/3 = 1.333..
d <- 20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))
set.seed(1)
n <- 200
U <- rnacopula(n, cop)
system.time(theta.hat.beta <- ebeta(U, cop=cop))
theta.hat.beta$root
system.time(theta.hat.dmle <- edmle(U, cop=cop))
theta.hat.dmle$minimum
system.time(theta.hat.etau <- etau(U, cop=cop, method="tau.mean"))
theta.hat.etau
system.time(theta.hat.etau. <- etau(U, cop=cop, method="theta.mean"))
theta.hat.etau.
## etau() in the case of missing values (NA's)
## ------ ---------------------
##' @title add Missing Values completely at random
##' @param x matrix or vector
##' @param prob desired probability ("fraction") of missing values (\code{\link{NA}}s).
##' @return x[] with some (100*prob percent) entries replaced by \code{\link{NA}}s.
addNAs <- function(x, prob) {
np <- length(x)
x[sample.int(np, prob*np)] <- NA
x
}
## UM[] := U[] with 5% missing values
set.seed(7); UM <- addNAs(U, p = 0.05)
mean(is.na(UM)) # 0.05
## This error if x has NA's was happening for etau(UM, cop=cop)
## before copula version 0.999-17 (June 2017) :
try(eM <- etau(UM, cop=cop, use = "everything"))
# --> Error ... NA/NaN/Inf in foreign function call
## The new default:
eM0 <- etau(UM, cop=cop, use = "pairwise")
eM1 <- etau(UM, cop=cop, use = "complete")
## use = "complete" is really equivalent to dropping all obs. with with missing values:
stopifnot(all.equal(eM1, etau(na.omit(UM), cop=cop), tol = 1e-15))
## but use = "pairwise" ---> cor(*, use = "pairwise") is much better:
rbind(etau.U = theta.hat.etau, etau.UM.pairwise = eM0, etau.UM.complete = eM1)