https://github.com/cran/nacopula
Raw File
Tip revision: 161411bb86f97e5a8bd89091cd61d03a33c2761a authored by Martin Maechler on 06 February 2012, 00:00:00 UTC
version 0.8-0
Tip revision: 161411b
cop_objects.R
## 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 <http://www.gnu.org/licenses/>.

#### 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)

back to top