## Copyright (C) 2012, 2015 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see .
###__ Estimation methods for Archimedean Copulas __
require(copula)
(doExtras <- copula:::doExtras())
## From source(system.file("test-tools.R", package = "Matrix")) :
showProc.time <- local({
pct <- proc.time()
function() { ## CPU elapsed __since last called__
ot <- pct ; pct <<- proc.time()
cat('Time elapsed: ', (pct - ot)[1:3],'\n')
}
})
### We have other implicit estimation "tests" in demo(estimation.gof)
## i.e., ../demo/estimation.gof.R
## ~~~~~~~~~~~~~~~~~~~~~~~~
## Frank & mle ---- log-density is *REALLY* not good enough:
## -----------
## at least not for "GoF" where we may have "tail-dependent"
## observations which are highly improbable under Frank and would
## "need" a very large theta :
p.log.f <- function(t0, t1, n.th = 256,
thet.vec = seq(t0, t1, length.out= n.th),
d, u0, uu = rbind(rep(u0, length.out = d)),
cop = "Frank", legend.x = "bottomright",
col = mapply(adjustcolor, col=c(1,3,2,4,1,2), alpha=c(2, 5:7, 4,2)/10),
lwd = c(5,4,3,1,5,5), lty = 1:6)
{
stopifnot(d == as.integer(d), length(d) == 1, d >= 2,
NCOL(uu) == d, is.numeric(thet.vec), length(thet.vec) > 1)
cop <- onacopulaL(getAcop(cop), list(NA, 1:d))
stopifnot(is.function(.f. <- cop@copula@dacopula))
## meths <- c("negI-s-Stirling", "negI-s-Eulerian")
meths <- local({m <- eval(formals(polylog)$method);m[grep("^negI-s", m)]})
l.arr <- sapply(meths, function(meth)
{
sapply(c(FALSE,TRUE), function(Li.log)
sapply(thet.vec, .f., u = uu, log=TRUE,
method = meth, Li.log.arg=Li.log))
}, simplify = "array")
dim(l.arr) <- local({D <- dim(l.arr); c(D[1], prod(D[-1]))})
## n x 4
nms <- c(t(outer(substring(meths, 8),
paste("Li_n(", c("", "log(arg)"), ")", sep=""),
##paste("Li.log.A=",c(FALSE,TRUE),sep=""),
paste, sep="_")))
colnames(l.arr) <- nms
tit <- bquote(.(cop@copula@name)*" copula - log density " ~~
log~f(theta, bold(u)[d==.(d)]))
matplot(thet.vec, l.arr, type = "l", lwd=lwd, col=col, lty=lty,
main = tit, xlab=expression(theta), ylab = expression(log~f(theta, .)))
legend(legend.x, nms, col=col, lwd=lwd, lty=lty,
title = "polylog(.., method = \"*\")", inset = .01)
mkU <- function(u, n1 = 4, n2 = 2) {
n <- length(u <- c(u))
cu <- if(n > (n1+n2)) ## use "..."
c(format(head(u, n1)), "...", format(tail(u, n2))) else format(u)
paste("(", paste(cu, collapse = ", "), ")", sep="")
}
mtext(bquote(bold(u)[d==.(d)] == .(mkU(uu))), line=1/8)
invisible(structure(cbind(theta = thet.vec, l.arr),
u = uu, dacopula = .f.))
}
unattr <- function(obj) { mostattributes(obj) <- NULL ; obj }
show. <- function(pmat) unattr(pmat)[, c(1,3,5,7)]
if(!dev.interactive(orNone=TRUE)) pdf("estim-ex.pdf")
th <- seq(1,60, by=1/4)
show.(m1 <- p.log.f(thet.vec = th, d = 5, u0 = .987))
## wow! asymptotic is quite good even here!
## using MC (instead of polylog(.)) also ``works'' {the function is noisy ..}:
l.th.MC <- sapply(th, attr(m1, "dacopula"), u = attr(m1, "u"),
n.MC = 2000, log=TRUE) ## 2000: be speedy in test
lines(th, l.th.MC, lwd=3, col=adjustcolor("sandybrown", .5))
legend("right", "n.MC = 2000", lwd=3, col=adjustcolor("sandybrown", .5),
inset=.01)
rm(th)
showProc.time()
## Extend the range:
show.(m2 <- p.log.f(1, 200, n.th=401, d = 5, u0 = .987))
if(doExtras) {
## Extend the range even more -- change u0
show.(m3 <- p.log.f(10, 500, d = 5, u0 = .96))
## and more
show.(m3.2 <- p.log.f(10, 800, d = 5, u0 = .96))#-> breakdown at ~ 775..
}
showProc.time()
## higher d:
if(doExtras) {
show.(m4 <- p.log.f(10, 500, d = 12, u0 = 0.95))
m5 <- p.log.f(10, 500, d = 12, u0 = 0.92)
m6 <- p.log.f(1, 200, d = 12, u0 = c(0.8, 0.9), legend.x="topright")
m7 <- p.log.f(1, 400, d = 12, u0 = c(0.88, 0.9))
}
## whereas this now also overflows for Eulerian *AND* asymp. + log.arg:
show.(mm <- p.log.f(10, 1000, d = 12, u0 = 0.9))
##--> investigation shows that w is *also* underflowing to 0 (!!!) :
myTracer <- quote({
cat(sprintf("Tracing %s: variables at entry are\n",
## the -4 is purely by trial and error -- better?
deparse(sys.call(-4)[1])))
print(ls.str())
})
trace(polylog, tracer = myTracer, print=FALSE)
## Tracing function "polylog" in package "copula"
f <- attr(mm,"dacopula")
f(attr(mm,"u"), 825, log=TRUE, method="negI-s-asymp")
## ....
## z : num -4.15e-322
## [1] 61.48259
f(attr(mm,"u"), 800, log=TRUE, method="negI-s-asymp")
## z : num -2.44e-312
if(FALSE)
debug(f)
## and then realize that in f(), i.e., Frank's dacopula() final line,
## (d-1)*log(theta) + Li. - theta*u.sum - lu
## the terms Li. - (theta * u.sum)
## 'more or less cancel' -- so if we could compute a
## *different* rescaled polylog() .. maybe we can get better
untrace(polylog)
if(doExtras)
## and similarly here:
ll <- p.log.f(1, 12000, d = 12, u0 = 0.08)
showProc.time()
###------- Diagonal MLE --- dDiag() , etc ----------------------------------
## Some basic dDiag() checks {that failed earlier}:
stopifnot(identical(0, dDiag(0, onacopulaL("AMH", list(0.5, 1:4)))),
all.equal(dDiag(.9, onacopulaL("Frank", list(44, 1:3))),
1.008252, tolerance= 1e-5),
all.equal(dDiag(.9, onacopulaL("Joe", list(44, 1:3))),
1.025283, tolerance= 1e-5),
TRUE)
demo("dDiag-plots", package = "copula")
## -----------> ../demo/dDiag-plots.R
showProc.time()
## A problem to solve:
curve(copAMH@dDiag(exp(-x), theta=.9, d=4, log=TRUE), 0, 500, n=2001,
col = 2, lwd= 2,
main = "FIXME (*not* urgent): dDiag(u, , log=TRUE) for small u")
showProc.time()
##----------- Test all estimation methods ----------------
## all methods
(estMeth <- eval(formals(enacopula)$method))
## Use selected (the best 6) estimation methods:
## (estMeth <- c("mle", "smle", "dmle", "tau.tau.mean", "mde.chisq.CvM", "mde.chisq.KS"))
n <- if(doExtras) 128 else 8
d <- 6 ; theta <- 2
(cop <- onacopula("Clayton", C(theta, 1:d)))
set.seed(1)
U <- rnacopula(n, cop)
## Show that 'racopula()' is indeed identical {and faster for the case of non-nested}
stopifnot(identical(U,
{set.seed(1); copula:::racopula(n, copClayton, theta=theta, d=d)}))
rr <- sapply(estMeth, function(e) {
enacopula(U, cop, e, n.MC = if(e == "smle") 1000 else 0)
})
round(cbind(rr, bias = rr - theta), 3)
showProc.time()
### 2-level nested Archimedean copulas :
if(doExtras) {
## was simply (!) : demo("dnac-demo")
v.dnac <- vignette("dNAC", package="copula")
R.dnac <-
if(nzchar(v.dnac$R))
with(v.dnac, file.path(Dir, "doc", R))
else { ## default: knitr::rmarkdown() unfortunately does *not* provide R
stopifnot(file.exists(ff <- with(v.dnac, file.path(Dir, "doc", File))))
oR <- file.path(tempdir(), sub("Rmd$", "R", basename(ff)))
knitr::purl(ff, output = oR)
oR
}
cat("R file: ", R.dnac, "\n -- now source()ing it :\n")
source(R.dnac, echo = TRUE, max.deparse.length = Inf,
keep.source = TRUE, encoding = getOption("encoding"))
}
showProc.time()