## Copyright (C) 2010--2011 Marius Hofert and Martin Maechler ## ## 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 . #### List of supported Archimedean copulas ## NOTA BENE: Write psi(), tau(), ... functions such that they *vectorize* ## --------- *and* do work even for (non-numeric) NULL argument ## {now checked in "acopula" validityMethod -- see ./AllClass.R } ### Ali-Mikhail-Haq, see Nelsen (2007) p. 116, # 3 ############################# copAMH <- (function() { ## to get an environment where C. itself is accessible C. <- new("acopula", name = "AMH", ## generator psi = function(t, theta) { (1-theta)/(exp(t+0)-theta) }, psiInv = function(u, theta, log=FALSE) { res <- log((1-theta*(1-u))/u) # alternative: log1p((1-theta)*(1/u-1)) if(log) log(res) else res }, ## parameter interval paraInterval = interval("[0,1)"), ## absolute value of generator derivatives psiDabs = function(t, theta, degree = 1, n.MC = 0, log = FALSE, is.log.t = FALSE, method = "negI-s-Eulerian", Li.log.arg=TRUE) { lth <- log(theta) if(n.MC > 0) { if(is.log.t) t <- exp(t) # very cheap for now psiDabsMC(t, family="AMH", theta=theta, degree=degree, n.MC=n.MC, log=log) } else { ## FIXME: deal with is.log.t ## Note: psiDabs(0, ...) is correct, namely (1-theta)/theta * polylog(theta, s=-degree) if(theta == 0) return(if(log) -t else exp(-t)) # independence Li.arg <- if(Li.log.arg) lth - t else theta*exp(-t) Li. <- polylog(Li.arg, s = -degree, method=method, is.log.z = Li.log.arg, log=log) if(log) Li. + log1p(-theta)-lth else Li. * (1-theta)/theta } }, ## derivatives of the generator inverse psiInvD1abs = function(u, theta, log=FALSE) { if(log) { log1p(-theta)-log(u)-log1p(-theta*(1-u)) } else { (1-theta)/(u*(1-theta*(1-u))) } }, ## density of the diagonal dDiag = function(u, theta, d, log=FALSE) { x <- (1-(1-u)*theta)/u if(any(iI <- is.infinite(x))) { ## for u == 0 (seems unneeded for small u << 1 ?) u[iI] <- if(log) -Inf else 0 ok <- !iI u[ok] <- C.@dDiag(u[ok], theta, d=d, log=log) return(u) } if(log) log(d)+ (d-1)*log(x) + 2*(log(x-theta)-log(x^d-theta)) else d* x^(d-1) * ((x-theta)/(x^d-theta))^2 }, ## density AMH dacopula = function(u, theta, n.MC=0, log=FALSE, method = "negI-s-Eulerian", Li.log.arg=TRUE) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 ## f() := NaN outside and on the boundary of the unit hypercube res <- rep.int(NaN, n <- nrow(u)) ## indices for which density has to be evaluated: n01 <- apply(u,1,function(x) all(0 < x, x < 1)) if(!any(n01)) return(res) if(theta == 0) { res[n01] <- if(log) 0 else 1; return(res) } # independence ## auxiliary results u. <- u[n01,, drop=FALSE] tIu <- -theta*(1-u.) sum. <- rowSums(log(u.)) sumIu <- rowSums(log1p(tIu)) ## main part if(n.MC > 0) { # Monte Carlo V <- C.@V0(n.MC, theta) res[n01] <- colMeans(exp(d*(log1p(-theta) + log(V)) + (V-1) %*% t(sum.) - (V+1) %*% t(sumIu))) if(log) log(res) else res } else { # explicit Li.arg <- if(Li.log.arg) log(theta) + sum. - sumIu else theta* apply(u./(1+tIu), 1, prod) Li. <- polylog(Li.arg, s = -d, method=method, is.log.z = Li.log.arg, log=TRUE) res[n01] <- (d+1)*log1p(-theta)-log(theta)+ Li. -sum. -sumIu if(log) res else exp(res) } }, ## score function score = function(u, theta) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 omu <- 1-u b <- rowSums(omu/(1-theta*omu)) h <- theta*apply(u/(1-theta*omu), 1, prod) -(d+1)/(1-theta) - 1/theta + b + (b+1/theta) * polylog(h, s=-(d+1), method="negI-s-Stirling") / polylog(h, s=-d, method="negI-s-Stirling") }, ## nesting constraint nestConstr = function(theta0,theta1) { C.@paraConstr(theta0) && C.@paraConstr(theta1) && theta1 >= theta0 }, ## V0 with density dV0 and V01 with density dV01 corresponding to ## LS^{-1}[exp(-V_0psi_0^{-1}(psi_1(t)))] V0 = function(n,theta) rgeom(n, 1-theta) + 1, dV0 = function(x,theta,log = FALSE) dgeom(x-1, 1-theta, log), V01 = function(V0,theta0,theta1) { rnbinom(length(V0),V0,(1-theta1)/(1-theta0))+V0 }, dV01 = function(x,V0,theta0,theta1,log = FALSE) { stopifnot(length(V0) == 1 || length(x) == length(V0)) dnbinom(x-V0,V0,(1-theta1)/(1-theta0),log = log) }, ## Kendall's tau tau = tauAMH, ##-> ./aux-acopula.R ## function(th) 1 - 2*((1-th)*(1-th)*log(1-th)+th)/(3*th*th) ## but numerically stable, including, theta -> 0 tauInv = function(tau, tol=.Machine$double.eps^0.25, check=TRUE, warn=TRUE, ...) { if((L <- any(tau > 1/3)) || any(tau < 0)) { ct <- if(length(ct <- sort(tau, decreasing = L)) > 3) paste(paste(format(ct[1:3]),collapse=", "),", ...",sep="") else format(ct) msg <- "For AMH, Kendall's tau must be in [0, 1/3], but (" if(check) stop(msg, if(L) "largest" else "smallest", " sorted) tau = ", ct) else { if(warn) warning(msg, if(L) "largest" else "smallest", " sorted) tau = ", ct) r <- tau r[sml <- tau <= 0 ] <- 0 r[lrg <- tau >= 1/3] <- 1 ok <- !sml & !lrg r[ok] <- C.@tauInv(tau[ok], tol=tol, ...) return(r) } } vapply(tau, function(tau) { if(abs(tau - 1/3) < 1e-10) return(1.) ## else: r <- safeUroot(function(th) tauAMH(th) - tau, interval = c(0, 1-1e-12), Sig = +1, tol = tol, check.conv=TRUE, ...) r$root }, 0.) }, ## lower tail dependence coefficient lambda_l lambdaL = function(theta) { 0*theta }, lambdaLInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for an Ali-Mikhail-Haq copula gives lambdaL = 0") NA * lambda }, ## upper tail dependence coefficient lambda_u lambdaU = function(theta) { 0*theta }, lambdaUInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for an Ali-Mikhail-Haq copula gives lambdaU = 0") NA * lambda } ) C. })()# {copAMH} ### Clayton, see Nelsen (2007) p. 116, #1 (slightly simpler form) ############## copClayton <- (function() { ## to get an environment where C. itself is accessible C. <- new("acopula", name = "Clayton", ## generator psi = function(t, theta) { (1+t)^(-1/theta) }, psiInv = function(u, theta, log=FALSE) { res <- u^(-theta) - 1 if(log) log(res) else res }, ## parameter interval paraInterval = interval("(0,Inf)"), ## absolute value of generator derivatives psiDabs = function(t, theta, degree=1, n.MC=0, log=FALSE) { if(n.MC > 0) { psiDabsMC(t, family="Clayton", theta=theta, degree=degree, n.MC=n.MC, log=log) } else { ## Note: psiDabs(0, ...) is correct, namely gamma(d+1/theta)/gamma(1/theta) alpha <- 1/theta res <- lgamma(degree+alpha)-lgamma(alpha)-(degree+alpha)*log1p(t) if(log) res else exp(res) } }, ## derivatives of the generator inverse psiInvD1abs = function(u, theta, log = FALSE) { if(log) log(theta)-(1+theta)*log(u) else theta*u^(-(1+theta)) }, ## density of the diagonal dDiag = function(u, theta, d, log=FALSE){ if(log) log(d)-(1+1/theta)*log(1+(d-1)*(1-u^theta)) else d*(1+(d-1)*(1-u^theta))^(-(1+1/theta)) }, ## density Clayton dacopula = function(u, theta, n.MC=0, log=FALSE) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 ## f() := NaN outside and on the boundary of the unit hypercube res <- rep.int(NaN, n <- nrow(u)) n01 <- apply(u,1,function(x) all(0 < x, x < 1)) # indices for which density has to be evaluated if(!any(n01)) return(res) ## auxiliary results u. <- u[n01,, drop=FALSE] lu <- rowSums(log(u.)) t <- rowSums(C.@psiInv(u., theta)) ## main part if(n.MC > 0) { # Monte Carlo lu.mat <- matrix(rep(lu, n.MC), nrow=n.MC, byrow=TRUE) V <- C.@V0(n.MC, theta) ## stably compute log(colMeans(exp(lx))) lx <- d*(log(theta) + log(V)) - log(n.MC) - (1+theta)*lu.mat - V %*% t(t) # matrix of exponents; dimension n.MC x n ["V x u"] ## note: smle goes wrong if: ## (1) d*log(theta*V) [old code] ## (2) U is small (simultaneously) ## (3) theta is large res[n01] <- lsum(lx) } else { # explicit res[n01] <- sum(log1p(theta*(0:(d-1)))) - (1+theta)*lu - (d+1/theta)*log1p(t) } if(log) res else exp(res) }, ## score function score = function(u, theta) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 lu <- log(u) t <- rowSums(C.@psiInv(u, theta=theta)) ltp1 <- log(1+t) lx <- t(log(-lu)-theta*lu) # caution: lsum needs an (d,n)-matrix ldt <- lsum(lx) # log of the derivative of t w.r.t. theta alpha <- 1/theta k <- 0:(d-1) sum(k/(theta*k+1))-rowSums(lu)+alpha^2*ltp1-(d+alpha)* exp(ldt-ltp1) }, ## nesting constraint nestConstr = function(theta0,theta1) { C.@paraConstr(theta0) && C.@paraConstr(theta1) && theta1 >= theta0 }, ## V0 with density dV0 and V01 with density dV01 corresponding to ## LS^{-1}[exp(-V_0psi_0^{-1}(psi_1(t)))] V0 = function(n,theta) { rgamma(n, shape = 1/theta) }, dV0 = function(x,theta,log = FALSE) dgamma(x, shape = 1/theta, log), V01 = function(V0,theta0,theta1) { retstable(alpha=theta0/theta1, V0) }, dV01 = function(x,V0,theta0,theta1,log = FALSE) { stopifnot(length(V0) == 1 || length(x) == length(V0)) alpha <- theta0/theta1 gamma <- (cos(pi/2*alpha)*V0)^(1/alpha) delta <- V0*(alpha == 1) ## NB: new dstable() is vectorized in (x, gamma, delta) [but not the others] dst <- dstable(x, alpha=alpha, beta = 1, gamma=gamma, delta=delta, pm = 1, log=log, tol = 128* .Machine$double.eps) if(log) V0-x + dst else exp(V0-x) * dst }, ## Kendall's tau tau = function(theta) { theta/(theta+2) }, tauInv = function(tau) { 2*tau/(1-tau) }, ## lower tail dependence coefficient lambda_l lambdaL = function(theta) { 2^(-1/theta) }, lambdaLInv = function(lambda) { -1/log2(lambda) }, ## upper tail dependence coefficient lambda_u lambdaU = function(theta) { 0*theta }, lambdaUInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for a Clayton copula gives lambdaU = 0") NA * lambda } ) C. })()# {copClayton} ### Frank, see Nelsen (2007) p. 116, # 5 ####################################### ##' Frank object copFrank <- (function() { ## to get an environment where C. itself is accessible C. <- new("acopula", name = "Frank", ## generator psi = function(t, theta) { -log1p(expm1(-theta)*exp(0-t))/theta ## == -log(1-(1-exp(-theta))*exp(-t))/theta }, psiInv = function(u, theta, log=FALSE) { ## == -log( (exp(-theta*u)-1) / (exp(-theta)-1) ) thu <- u*theta # (-> recycling args) if(!length(thu)) return(thu) # {just for numeric(0) ..hmm} et1 <- expm1(-theta) # e^{-th} - 1 < 0 ## FIXME ifelse() is not quite efficient ## FIXME(2): the "> c* th" is pi*Handgelenk ### FIXME: use delta = exp(-thu)*(1 - exp(thu-th)/ (-et1) = ### = exp(-thu)* expm1(thu-theta)/et1 (4) ###-- do small Rmpfr-study to find the best form -- (4) above and the three forms below r <- ifelse(thu > .01*theta, # thu = u*th > .01*th <==> u > 0.01 { e.t <- exp(-theta) ifelse(e.t > 0 & abs(theta - thu) < 1/2,# th -th*u = th(1-u) < 1/2 -log1p(e.t * expm1(theta - thu)/et1), -log1p((exp(-thu)- e.t) / et1)) }, ## for small u (u < 0.01) : -log(expm1(-thu)/et1)) if(log) log(r) else r }, ## parameter interval paraInterval = interval("(0,Inf)"), ## absolute value of generator derivatives psiDabs = function(t, theta, degree = 1, n.MC = 0, log = FALSE, is.log.t = FALSE, method = "negI-s-Eulerian", Li.log.arg = TRUE) { if(n.MC > 0) { psiDabsMC(t, family="Frank", theta=theta, degree=degree, n.MC=n.MC, log=log) } else { ## Note: psiDabs(0, ...) is correct, namely (1/theta)*polylog(1-exp(-theta), s=-(degree-1)) Li.arg <- if(Li.log.arg) log1mexpm(theta) - t else -expm1(-theta)*exp(-t) Li. <- polylog(Li.arg, s = -(degree-1), log=log, method=method, is.log.z = Li.log.arg) if(log) Li. - log(theta) else Li. / theta } }, ## derivatives of the generator inverse psiInvD1abs = function(u, theta, log = FALSE) { if(log) log(theta)- {y <- u*theta; y + log1mexpm(y)} else theta/expm1(u*theta) }, ## density of the diagonal dDiag = function(u, theta, d, log=FALSE) { r <- ut <- u*theta ## h <- -expm1(-ut) # 1 - exp(-u * theta) ## L <- ut > log(2)*.Machine$double.digits # <==> exp(-ut) < Mach..eps <==> h=1 ## empirically, this is "better:" ## (more in ../demo/dDiag-plots.R -- should also depend on d ! L <- ut > 25 # [L]arge S <- ut < 0.10 # [S]mall M <- !L & !S # [M]edium ## recycle (u, theta): u <- rep(u, length.out=length(ut)) th <- rep(theta, length.out=length(ut)) r[L] <- dDiagFrank(u[L], th[L], d, log=log, method = "poly2") r[M] <- dDiagFrank(u[M], th[M], d, log=log, method = "poly4") r[S] <- dDiagFrank(u[S], th[S], d, log=log, method = "m1") r }, ## density Frank dacopula = function(u, theta, n.MC=0, log=FALSE, method = "negI-s-Eulerian", Li.log.arg=TRUE) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 ## f() := NaN outside and on the boundary of the unit hypercube res <- rep.int(NaN, n <- nrow(u)) n01 <- apply(u,1,function(x) all(0 < x, x < 1)) # indices for which density has to be evaluated if(!any(n01)) return(res) ## auxiliary results u. <- u[n01,, drop=FALSE] u.sum <- rowSums(u.) lp <- log1mexpm(theta) # log(1 - exp(-theta)) lpu <- log1mexpm(theta*u.) # log(1 - exp(-theta * u)) lu <- rowSums(lpu) ## main part res[n01] <- if(n.MC > 0) { # Monte Carlo V <- C.@V0(n.MC, theta) lx <- rep(-theta*u.sum, each=n.MC) + d*(log(theta) + log(V) - V*lp) - log(n.MC) + (V-1) %*% t(lu) lsum(lx) } else { # explicit Li.arg <- if(Li.log.arg) lp + rowSums(lpu-lp) else -expm1(-theta)*exp(rowSums(lpu-lp)) Li. <- polylog(Li.arg, s = -(d-1), log=TRUE, method=method, is.log.z = Li.log.arg) (d-1)*log(theta) + Li. - theta*u.sum - lu } if(log) res else exp(res) }, ## score function score = function(u, theta) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 e <- exp(-theta) Ie <- -expm1(-theta) # == 1 - e == 1-e^{-theta} etu <- exp(mtu <- -theta*u) # exp(-theta*u) Ietu <- -expm1(mtu) ## = 1 - etu == 1 - exp(-theta*u) ## FIXME: allow 'Li.log.arg' -> psilog(*, is.log.z) as for dacopula() above h <- Ie*apply(Ie/etu, 1, prod) factor <- rowSums(u*etu/Ietu) - (d-1)*e/Ie (d-1)/theta - rowSums(u/Ietu) + factor * polylog(h, s=-d, method="negI-s-Stirling") / polylog(h, s=-(d-1), method="negI-s-Stirling") }, ## nesting constraint nestConstr = function(theta0,theta1) { C.@paraConstr(theta0) && C.@paraConstr(theta1) && theta1 >= theta0 }, ## V0 (algorithm of Kemp (1981)) with density dV0 and V01 with density ## dV01 corresponding to LS^{-1}[exp(-V_0psi_0^{-1}(psi_1(t)))] V0 = function(n,theta) rlog(n, -expm1(-theta), exp(-theta)), dV0 = function(x,theta, log = FALSE) { if(any(x != (x. <- round(x)))) { x <- x.; warning("x has been rounded to integer") } if(log) { x*log1mexpm(theta) - log(x*theta) } else { p <- -expm1(-theta) # == 1-exp(-theta) p^x/(x*theta) } }, V01 = function(V0, theta0, theta1, rej = 1, approx = 10000) { ## rej method switch: if V0*theta_0*p0^(V0-1) > rej a rejection ## from F01 of Joe is applied (otherwise, the sum is ## sampled via a logarithmic envelope for the summands) ## approx is the largest number of summands before asymptotics is used (see copJoe@V01) ## FIXME: optimal value of rej (for approx = 10000) is not clear yet; rej = 1 is not bad, however; ## lgammacor gives underflow warnings for rej < 1 rF01Frank(V0, theta0, theta1, rej, approx) }, dV01 = function(x,V0,theta0,theta1,log = FALSE) { stopifnot(length(V0) == 1 || length(x) == length(V0), all(x >= V0)) lfactor <- x*log1mexpm(theta1) - V0*log1mexpm(theta0) res <- lfactor + dsumSibuya(x, V0, theta0/theta1, log=TRUE) if(log) res else exp(res) }, ## Kendall's tau; debye_1() is from package 'gsl' : tau = function(theta) { if((l <- length(theta)) == 0) return(numeric(0)) # to work with NULL res <- numeric(l) res[isN <- theta == 0] <- 0 # limiting case res[na <- is.na(theta)] <- NA if(any(i <- !(na | isN))) res[i] <- 1 + 4*(debye_1(theta[i]) - 1)/theta[i] res }, tauInv = function(tau, tol = .Machine$double.eps^0.25, ...) { res <- tau res[isN <- res == 0] <- 0 # limiting case res[!isN] <- sapply(res[!isN], function(tau) { r <- safeUroot(function(th) C.@tau(th) - tau, interval = c(1e-12, 198), Sig = +1, tol=tol, check.conv=TRUE, ...) r$root }) res }, ## lower tail dependence coefficient lambda_l lambdaL = function(theta) { 0*theta }, lambdaLInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for a Frank copula gives lambdaL = 0") NA * lambda }, ## upper tail dependence coefficient lambda_u lambdaU = function(theta) { 0*theta }, lambdaUInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for a Frank copula gives lambdaU = 0") NA * lambda } ) C. })()# {copFrank} ### Gumbel, see Nelsen (2007) p. 116, # 4 ###################################### copGumbel <- (function() { ## to get an environment where C. itself is accessible C. <- new("acopula", name = "Gumbel", ## generator psi = function(t, theta) { exp(-t^(1/theta)) }, psiInv = function(u, theta, log=FALSE) { if(log) theta*log(-log(u)) else (-log(u+0))^theta }, ## parameter interval paraInterval = interval("[1,Inf)"), ## absolute value of generator derivatives psiDabs = function(t, theta, degree=1, n.MC=0, method = eval(formals(polyG)$method), log = FALSE) { is0 <- t == 0 isInf <- is.infinite(t) res <- numeric(n <- length(t)) res[is0] <- if(theta==1) 0 else Inf # Note: psiDabs(0, ...) is correct (even for n.MC > 0) res[isInf] <- -Inf n0Inf <- !(is0 | isInf) if(all(!n0Inf)) return(if(log) res else exp(res)) t. <- t[n0Inf] if(n.MC > 0) { res[n0Inf] <- psiDabsMC(t, family="Gumbel", theta=theta, degree=degree, n.MC=n.MC, log=TRUE) } else { if(theta == 1) { res[n0Inf] <- -t. # independence } else { alpha <- 1/theta lt <- log(t.) res[n0Inf] <- -degree*lt -t.^alpha + polyG(alpha*lt, alpha=alpha, d = degree, method=method, log=TRUE) } } if(log) res else exp(res) }, ## derivatives of the generator inverse psiInvD1abs = function(u, theta, log = FALSE) { lu <- log(u) if(log) { log(theta)+(theta-1)*log(-lu)-lu } else { theta*(-lu)^(theta-1)/u } }, ## density of the diagonal dDiag = function(u, theta, d, log=FALSE) { alpha <- 1/theta da <- d^alpha if(log) (da-1)*log(u) + alpha*log(d) else da*u^(da-1) }, ## density Gumbel dacopula = function(u, theta, n.MC=0, method = eval(formals(polyG)$method), log=FALSE) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 ## f() := NaN outside and on the boundary of the unit hypercube res <- rep.int(NaN, n <- nrow(u)) n01 <- apply(u,1,function(x) all(0 < x, x < 1)) # indices for which density has to be evaluated if(!any(n01)) return(res) if(theta == 1) { res[n01] <- if(log) 0 else 1; return(res) } # independence ## auxiliary results u. <- u[n01,, drop=FALSE] mlu <- -log(u.) # -log(u) lmlu <- log(mlu) # log(-log(u)) t <- rowSums(C.@psiInv(u., theta)) ## main part if(n.MC > 0) { # Monte Carlo V <- C.@V0(n.MC, theta) sum. <- rowSums((theta-1)*lmlu + mlu) sum.mat <- matrix(rep(sum., n.MC), nrow=n.MC, byrow=TRUE) ## stably compute log(colMeans(exp(lx))) lx <- d*(log(theta)+log(V)) - log(n.MC) - V %*% t(t) + sum.mat # matrix of exponents; dimension n.MC x n ["V x u"] res[n01] <- lsum(lx) if(log) res else exp(res) } else { # explicit alpha <- 1/theta ## compute lx = alpha*log(sum(psiInv(u., theta))) lx <- alpha*log(t) ## ==== former version [start] (numerically slightly more stable but slower) ==== ## im <- apply(u., 1, which.max) ## mat.ind <- cbind(seq_len(n), im) # indices that pick out maxima from u. ## mlum <- mlu[mat.ind] # -log(u_max) ## mlum.mat <- matrix(rep(mlum, d), ncol = d) ## lx <- lmlu[mat.ind] + alpha*log(rowSums((mlu/mlum.mat)^theta)) # alpha*log(sum(psiInv(u, theta))) ## ==== former version [end] ==== ## compute sum ls. <- polyG(lx, alpha, d, method=method, log=TRUE)-d*lx/alpha ## the rest cop.val <- pnacopula(onacopulaL("Gumbel", list(theta, 1:d)), u.) res[n01] <- d*log(theta) + rowSums((theta-1)*lmlu + mlu) + ls. res[n01] <- if(log) log(cop.val) + res[n01] else cop.val * exp(res[n01]) res } }, ## score function score = function(u, theta) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 stop("The score function is currently not implemented for Gumbel copulas") }, ## nesting constraint nestConstr = function(theta0,theta1) { C.@paraConstr(theta0) && C.@paraConstr(theta1) && theta1 >= theta0 }, ## V0 with density dV0 and V01 with density dV01 corresponding to ## LS^{-1}[exp(-V_0psi_0^{-1}(psi_1(t)))] V0 = function(n,theta) { if(theta == 1) { ## Sample from S(1,1,0,1;1) ## with Laplace-Stieltjes transform exp(-t) rep.int(1., n) } else { alpha <- 1/theta ## Sample from S(alpha,1,(cos(alpha*pi/2))^(1/alpha),0;1) ## with Laplace-Stieltjes transform exp(-t^alpha) rstable1(n, alpha, beta=1, gamma = cos(alpha*pi/2)^(1/alpha)) ## Note: calling sequence: ## rstable1 == rstable1C (in rstable1.R) ## -> rstable_c (in retstable.c) -> rstable_vec (in retstable.c) ## -> rstable0 (in retstable.c) } }, dV0 = function(x,theta,log = FALSE) C.@dV01(x,1,1,theta,log), V01 = function(V0,theta0,theta1) { alpha <- theta0/theta1 if(alpha == 1) { ## Sample from S(1,1,0,V0;1) ## with Laplace-Stieltjes transform exp(-V0*t) V0 } else { rstable1(length(V0), alpha, beta=1, gamma = (cos(alpha*pi/2)*V0)^(1/alpha)) ## Sample from S(alpha,1,(cos(alpha*pi/2)V0)^(1/alpha),0;1) ## with Laplace-Stieltjes transform exp(-V0*t^alpha) } }, dV01 = function(x,V0,theta0,theta1,log = FALSE) { stopifnot(length(V0) == 1 || length(x) == length(V0)) alpha <- theta0/theta1 gamma <- (cos(pi/2*alpha)*V0)^(1/alpha) delta <- V0*(alpha == 1) ## NB: new dstable() is vectorized in (x, gamma, delta) [but not the others] dstable(x, alpha=alpha, beta = 1, gamma=gamma, delta=delta, pm = 1, log=log, tol = 128* .Machine$double.eps) }, ## Kendall's tau tau = function(theta) { (theta-1)/theta }, tauInv = function(tau) { 1/(1-tau) }, ## lower tail dependence coefficient lambda_l lambdaL = function(theta) { 0*theta }, lambdaLInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for a Gumbel copula gives lambdaL = 0") NA * lambda }, ## upper tail dependence coefficient lambda_u lambdaU = function(theta) { 2 - 2^(1/theta) }, lambdaUInv = function(lambda) { 1/log2(2-lambda) } ) C. })()# {copGumbel} ### Joe, see Nelsen (2007) p. 116, # 6 ######################################### ##' Joe object copJoe <- (function() { ## to get an environment where C. itself is accessible C. <- new("acopula", name = "Joe", ## generator psi = function(t, theta) { 1 - (-expm1(0-t))^(1/theta) ## == 1 - (1-exp(-t))^(1/theta) }, psiInv = function(u, theta, log=FALSE) { res <- -log1p(-(1-u)^theta) if(log) log(res) else res }, ## parameter interval paraInterval = interval("[1,Inf)"), ## absolute value of generator derivatives psiDabs = function(t, theta, degree = 1, n.MC = 0, method= eval(formals(polyJ)$method), log = FALSE) { is0 <- t == 0 isInf <- is.infinite(t) res <- numeric(n <- length(t)) res[is0] <- if(theta==1) 0 else Inf # Note: psiDabs(0, ...) is correct (even for n.MC > 0) res[isInf] <- -Inf n0Inf <- !(is0 | isInf) if(all(!n0Inf)) return(if(log) res else exp(res)) t. <- t[n0Inf] if(n.MC > 0) { res[n0Inf] <- psiDabsMC(t, family="Joe", theta=theta, degree=degree, n.MC=n.MC, log=TRUE) } else { if(theta == 1) { res[n0Inf] <- -t. # independence } else { alpha <- 1/theta mt <- -t. l1mt <- log1mexpm(t.) # log(1-exp(-t)) sum. <- polyJ(mt-l1mt, alpha, degree, method=method, log=TRUE) res[n0Inf] <- -log(theta) + mt - (1-alpha)*l1mt + sum. } } if(log) res else exp(res) }, ## derivatives of the generator inverse psiInvD1abs = function(u, theta, log = FALSE) { if(log) { log(theta)+(theta-1)*log1p(-u)-log1p(-(1-u)^theta) } else { theta/((1-u)^(1-theta)-(1-u)) } }, ## density of the diagonal dDiag = function(u, theta, d, log=FALSE) { ## d* x^(d-1) * ((1-x^d)/(1-x)) ^ (1/theta-1), where x = 1 - (1-u)^theta ## we use (1-x^d)/(1-x) === circRat((1-u)^theta, d) Ix <- (1-u)^theta x <- 1-Ix if(log) log(d)+ (d-1)*log(x) + (1/theta-1)*log(circRat(Ix, d)) ## FIXME? for log-case: log(circRat(Ix, d)) = log1p(-x^d)-theta*log1p(-u)) else d* x^(d-1) * circRat(Ix, d)^(1/theta-1) }, ## density Joe dacopula = function(u, theta, n.MC=0, method = eval(formals(polyJ)$method), log = FALSE) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 ## f() := NaN outside and on the boundary of the unit hypercube res <- rep.int(NaN, n <- nrow(u)) n01 <- apply(u,1,function(x) all(0 < x, x < 1)) # indices for which density has to be evaluated if(!any(n01)) return(res) if(theta == 1) { res[n01] <- if(log) 0 else 1; return(res) } # independence ## auxiliary results u. <- u[n01,, drop=FALSE] l1_u <- rowSums(log1p(-u.)) # sum_j log(1-u_j) lh <- rowSums(log1p(-(1-u.)^theta)) # rowSums(log(1-(1-u)^theta)) = log(h) ## main part if(n.MC > 0) { # Monte Carlo V <- C.@V0(n.MC, theta) sum. <- (theta-1)*l1_u sum.mat <- matrix(rep(sum., n.MC), nrow=n.MC, byrow=TRUE) ## stably compute log(colMeans(exp(lx))) ## matrix of exponents; dimension n.MC x n ["V x u"] : lx <- d*(log(theta) + log(V)) + (V-1) %*% t(lh) + sum.mat - log(n.MC) res[n01] <- lsum(lx) } else { alpha <- 1/theta l1_h <- log1mexpm(-lh) # log(1-h) lh_l1_h <- lh - l1_h # log(h/(1-h)) res[n01] <- (d-1)*log(theta) + (theta-1)*l1_u - (1-alpha)*l1_h + polyJ(lh_l1_h, alpha, d, method=method, log=TRUE) } if(log) res else exp(res) }, ## score function score = function(u, theta, method=eval(formals(polyJ)$method)) { stopifnot(C.@paraConstr(theta)) if(!is.matrix(u)) u <- rbind(u) if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2 l1_u <- rowSums(log1p(-u)) # log(1-u) u.th <- (1-u)^theta # (1-u)^theta lh <- rowSums(log1p(-u.th)) # rowSums(log(1-(1-u)^theta)) = log(h) l1_h <- log1mexpm(-lh) # log(1-h) lh_l1_h <- lh - l1_h # log(h/(1-h)) b <- rowSums(-l1_u*u.th/(1-u.th)) lP <- polyJ(lh_l1_h, alpha, d, method=method, log=TRUE) k <- 1:d alpha <- 1/theta s <- alpha * unlist(lapply(k, function(k.) sum(1/(theta*(1:k.)-1)))) ls <- log(s. + (k-1)*exp(log(b) - l1_h)) l.a.k <- log(Stirling2.all(d)) + lgamma(k-alpha) - lgamma(1-alpha) + ls lQ <- lsum(l.a.k + (k-1) %*% t(lh_l1_h)) (d-1)/theta + rowSums(l1_u) - l1_h/theta^2 + (1-1/theta)* lh_l1_h*b + exp(lQ-lP) }, ## nesting constraint nestConstr = function(theta0,theta1) { C.@paraConstr(theta0) && C.@paraConstr(theta1) && theta1 >= theta0 }, ## V0 with density dV0 and V01 with density dV01 corresponding to ## LS^{-1}[exp(-V_0psi_0^{-1}(psi_1(t)))] V0 = function(n,theta) rSibuya(n, 1/theta), dV0 = function(x,theta,log = FALSE) { if(log) lchoose(1/theta,x) else abs(choose(1/theta,x)) }, V01 = function(V0, theta0, theta1, approx = 10000) { ## approx is the largest number of summands before asymptotics is used alpha <- theta0/theta1 rF01Joe(V0, alpha, approx) }, dV01 = function(x, V0, theta0, theta1, method= eval(formals(dsumSibuya)$method), log = FALSE) { stopifnot(length(V0) == 1 || length(x) == length(V0)) ## also holds for theta0 == theta1 ## note: this is numerically challenging dsumSibuya(x, V0, theta0/theta1, method=method, log=log) }, ## Kendall's tau ## noTerms: even for theta==0, the approximation error is < 10^(-5) tau = function(theta, noTerms=446) { k <- noTerms:1 sapply(theta, function(th) { tk2 <- th*k + 2 1 - 4*sum(1/(k*tk2*(tk2 - th))) ## ==... (1/(k*(th*k+2)*(th*(k-1)+2))) }) }, tauInv = function(tau, tol = .Machine$double.eps^0.25, ...) { sapply(tau,function(tau) { r <- safeUroot(function(th) C.@tau(th) - tau, interval = c(1, 98), Sig = +1, tol=tol, check.conv=TRUE, ...) r$root }) }, ## lower tail dependence coefficient lambda_l lambdaL = function(theta) { 0*theta }, lambdaLInv = function(lambda) { if(any(lambda != 0)) stop("Any parameter for a Joe copula gives lambdaL = 0") NA * lambda }, ## upper tail dependence coefficient lambda_u lambdaU = function(theta) { 2-2^(1/theta) }, lambdaUInv = function(lambda) { log(2)/log(2-lambda) } ) C. })()# {copJoe} ### naming stuff ############################################################### cNms <- c("copAMH", "copClayton", "copFrank", "copGumbel", "copJoe") ## == dput(ls("package:nacopula",pat="^cop")) nmsC <- unlist(lapply(cNms, function(.)get(.)@name)) sNms <- abbreviate(nmsC, 1) ## keep these {hidden, for now}: c_shortNames <- structure(sNms, names = nmsC) c_longNames <- structure(nmsC, names = sNms) c_objNames <- structure(cNms, names = nmsC) rm(cNms, nmsC, sNms)