## Copyright (C) 2010 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 . ##' Outer power transformation of Archimedean copulas ##' ##' @title Outer power Archimedean copulas ##' @param copbase a "base" copula, i.e. of class "acopula"; ##' must be one of the predefined families ##' @param thetabase the (univariate) parameter 'theta' for the base copula ##' @return a new "acopula" object; the outer power copula [with generator psi(t^(1/theta)] ##' @author Marius Hofert opower <- function(copbase, thetabase) { ## create object with name in here so it's in the environment and we can access it C. <- new("acopula", name = paste("opower", copbase@name, sep=":"), ## generator psi = function(t,theta) { copbase@psi(t^(1/theta), thetabase) }, psiInv = function(t,theta, log=FALSE) { if(log) theta * copbase@psiInv(t, thetabase, log=TRUE) else copbase@psiInv(t, thetabase)^theta }, ## parameter interval paraInterval = interval("[1,Inf)"), ## nesting constraint nestConstr = function(theta0,theta1) { copbase@paraConstr(theta0) && copbase@paraConstr(theta1) && theta1 >= theta0 }, ## absolute value of generator derivatives psiDabs = function(t, theta, degree=1, n.MC=0, method=c("stirling", "binomial.coeff"), log=FALSE) { if(theta == 1) return(copbase@psiDabs(t, theta, degree=degree, n.MC=n.MC, log=log)) # copbase case is0 <- t == 0 isInf <- is.infinite(t) res <- numeric(n <- length(t)) res[is0] <- 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=C., theta=theta, degree=degree, n.MC=n.MC, log=TRUE) } else { k <- 1:degree # compute everything once for 1:degree beta <- 1/theta t.beta <- t.^beta # beta = 1/theta for psi(t^beta) ## in principle, it would be efficient to vectorize ## the psiDabs slots also in the parameter "degree", ## but that would be even more complicated lpsiDabs <- do.call(rbind, lapply(k, function(k.) copbase@psiDabs(t.beta, theta=thetabase, degree=k., log=TRUE))) # (degree,n)-matrix bklt <- k %*% t(log(t.beta)) # (degree,n)-matrix method <- match.arg(method) res[n0Inf] <- switch(method, "stirling" = { # much faster & less prone to errors s <- Stirling1.all(degree) b.one.j <- function(j.) { k <- 1:j. signs <- (-1)^k lS <- log(Stirling2.all(j.)) a <- lS + lpsiDabs[k,, drop=FALSE] + bklt[k,, drop=FALSE] (-beta)^j. * abs(s[j.]) * colSums(signs*exp(a)) } ## => returns a vector of length n containing the values for one j. and all t j <- 1:degree b <- do.call(rbind, lapply(j, FUN=b.one.j)) # (degree, n)-matrix -degree*log(t.)+log(colSums(b)) }, "binomial.coeff" = { ## outer sum lfac <- lfactorial(0:degree) # log(0!), log(1!), .., log(degree!) log.b.one.j <- function(j.) { k <- j.:degree a <- lpsiDabs[k,, drop=FALSE] + bklt[k,, drop=FALSE] - lfac[j.+1] - lfac[k-j.+1] # (degree-j.+1, n)-matrix ls. <- lsum(a) # length = n lchoose(beta*j., degree) + ls. # note: the lchoose() can be non-finite with this approach! } ## => returns a vector of length n containing the values for one j. and all t j <- 1:degree b <- do.call(rbind, lapply(j, FUN=log.b.one.j)) # (degree, n)-matrix signs <- signFF(beta, j, degree) lfac[degree+1] - degree*log(t.) + lssum(b, signs, strict=FALSE) }, stop(sprintf("unsupported method '%s' in psiDabs", method))) # end{switch} } if(log) res else exp(res) }, ## derivatives of the generator inverse psiInvD1abs = function(t, theta, log=FALSE) { if(theta == 1) return(copbase@psiInvD1abs(t, theta, log=log)) # copbase case if(log) { log(theta)+(theta-1)*log(copbase@psiInv(t,thetabase))+ copbase@psiInvD1abs(t, thetabase,log=TRUE) } else { theta*copbase@psiInv(t,thetabase)^(theta-1)* copbase@psiInvD1abs(t, thetabase,log=FALSE) } }, ## density dacopula = function(u, theta, n.MC=0, log=FALSE) { stopifnot(C.@paraConstr(theta)) if(theta == 1) return(copbase@dacopula(u, theta, n.MC=n.MC, log=log)) # copbase case 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] psiI <- rowSums(C.@psiInv(u.,theta)) res[n01] <- C.@psiDabs(psiI, theta, degree=d, n.MC=n.MC, log=TRUE) + rowSums(C.@psiInvD1abs(u., theta, log=TRUE)) 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 stop("The score function is currently not implemented for outer power copulas") }, ## V0 and V01 V0 = function(n, theta) { V0base <- copbase@V0(n, thetabase) if(theta == 1) return(V0base) # the copula is copbase with thetabase alpha <- 1/theta ## Sample from S(alpha,1,(cos(alpha*pi/2))^(1/alpha),0;1) ## with Laplace-Stieltjes transform exp(-t^alpha) S <- rstable1(n, alpha, beta=1, gamma = (cos(alpha*pi/2))^(1/alpha)) S*V0base^theta }, dV0 = function(x, theta, log=FALSE) { stop("not implemented; it's the density of SV^theta, where V ~ F with LS[F] = copbase@psi and S ~ S(1/theta, 1, cos^theta(pi/(2*theta)), I_{theta==1}; 1)") }, 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) { copGumbel@dV01(x, V0, theta0, theta1, log=log) }, ## Kendall's tau tau = function(theta) { 1-(1-copbase@tau(thetabase))/theta }, tauInv = function(tau) { taubase <- copbase@tau(thetabase) if(tau >= taubase) (1-taubase)/(1-tau) else { stop("The provided tau has to be >= taubase") NA * tau } }, ## lower tail dependence coefficient lambda_l lambdaL = function(theta) { if(copbase@name=="Clayton") 2^(-1/(thetabase*theta)) else 0*theta }, lambdaLInv = function(lambda) { if(copbase@name=="Clayton") { if(lambda >= 2^(-1/thetabase)) -1/(thetabase*log2(lambda)) else { stop("The provided lambda has to be >= 2^(-1/thetabase)") NA * lambda } } else { if(any(lambda != 0)) stop("Any parameter for this outer power Archimedean copula gives lambdaL = 0") NA * lambda } }, ## upper tail dependence coefficient lambda_u lambdaU = function(theta) { 2 - 2^(1/if(copbase@name %in% c("Gumbel", "Joe")) thetabase*theta else theta) }, lambdaUInv = function(lambda) { if(copbase@name %in% c("Gumbel", "Joe")) { if(lambda >= 2-2^(1/thetabase)) 1/(thetabase*log2(2-lambda)) else { stop("The provided lambda has to be >= 2-2^(1/thetabase)") NA * lambda } } else 1/log2(2-lambda) } ) C. }