https://github.com/cran/nacopula
Revision 161411bb86f97e5a8bd89091cd61d03a33c2761a authored by Martin Maechler on 06 February 2012, 00:00:00 UTC, committed by Gabor Csardi on 06 February 2012, 00:00:00 UTC
1 parent 5bc804b
Tip revision: 161411b
trafos.R
``````## Copyright (C) 2010 Marius Hofert and Martin Maechler
##
## This program is free software; you can redistribute it and/or modify it under
## 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/>.

##' 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.
}
``````

Computing file changes ...