## 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 . #### List of supported Archimedean copulas ## FIXME: Not "nice" that the names of the copula objects must be be used ## inside some of their components. ## Possible solution: ## use *same* environment for (tau, tauInv, paraConstr, nestConstr) ## 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 <- new("acopula", name = "AMH", ## generator psi = function(t,theta) { (1-theta)/(exp(t+0)-theta) }, psiInv = function(t,theta) { log((1-theta*(1-t))/t) }, ## parameter interval paraInterval = interval("[0,1)"), ## nesting constraint nestConstr = function(theta0,theta1) { copAMH@paraConstr(theta0) && copAMH@paraConstr(theta1) && theta1 >= theta0 }, ## V0 and V01 V0 = function(n,theta) rgeom(n, 1-theta) + 1, V01 = function(V0,theta0,theta1) { rnbinom(length(V0),V0,(1-theta1)/(1-theta0))+V0 }, ## 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, ...) { if(any(tau > 1/3)) stop("Impossible for AMH copula to attain a Kendall's tau larger than 1/3") sapply(tau,function(tau) { r <- safeUroot(function(th) tauAMH(th) - tau, interval = c(1e-12, 1-1e-12), 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 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 }) ### ==== Clayton, see Nelsen (2007) p. 116, #1 (slightly simpler form here) ==== copClayton <- new("acopula", name = "Clayton", ## generator psi = function(t,theta) { (1+t)^(-1/theta) }, psiInv = function(t,theta) { t^(-theta) - 1 }, ## parameter interval paraInterval = interval("(0,Inf)"), ## nesting constraint nestConstr = function(theta0,theta1) { copClayton@paraConstr(theta0) && copClayton@paraConstr(theta1) && theta1 >= theta0 }, ## V0 and V01 V0 = function(n,theta) { rgamma(n, shape = 1/theta) }, V01 = function(V0,theta0,theta1) { retstable(alpha=theta0/theta1, V0) }, ## 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 }) ### ==== Frank, see Nelsen (2007) p. 116, # 5 ================================== ##' Frank object copFrank <- 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(t,theta) { -log(expm1(-theta*t)/expm1(-theta)) ## == -log((exp(-theta*t)-1)/(exp(-theta)-1)) }, ## parameter interval paraInterval = interval("(0,Inf)"), ## nesting constraint nestConstr = function(theta0,theta1) { copFrank@paraConstr(theta0) && copFrank@paraConstr(theta1) && theta1 >= theta0 }, ## V0 (algorithm of Kemp (1981)) and V01 V0 = function(n,theta) { rlog(n,1-exp(-theta)) }, V01 = function(V0,theta0,theta1) { ## FIXME: how to approximate when V0 large? not yet solved ## theoretically sapply(lapply(V0, rFFrank, theta0=theta0,theta1=theta1), sum) }, ## Kendall's tau; debye_1() is from package 'gsl' : tau = function(theta) 1 + 4*(debye_1(theta) - 1)/theta, tauInv = function(tau, tol = .Machine$double.eps^0.25, ...) { sapply(tau, function(tau) { r <- safeUroot(function(th) copFrank@tau(th) - tau, interval = c(0.001,100), 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 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 }) ### ==== Gumbel, see Nelsen (2007) p. 116, # 4 ================================= copGumbel <- new("acopula", name = "Gumbel", ## generator psi = function(t,theta) { exp(-t^(1/theta)) }, psiInv = function(t,theta) { (-log(t+0))^theta }, ## parameter interval paraInterval = interval("[1,Inf)"), ## nesting constraint nestConstr = function(theta0,theta1) { copGumbel@paraConstr(theta0) && copGumbel@paraConstr(theta1) && theta1 >= theta0 }, ## V0 and V01 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)) } }, 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) } }, ## 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) } ) # {copGumbel} ### ==== Joe, see Nelsen (2007) p. 116, # 6 ==================================== ##' Joe object copJoe <- new("acopula", name = "Joe", ## generator psi = function(t,theta) { 1 - (-expm1(0-t))^(1/theta) ## == 1 - (1-exp(-t))^(1/theta) }, psiInv = function(t,theta) { -log1p(-(1-t)^theta) }, ## parameter interval paraInterval = interval("[1,Inf)"), ## nesting constraint nestConstr = function(theta0,theta1) { copJoe@paraConstr(theta0) && copJoe@paraConstr(theta1) && theta1 >= theta0 }, ## V0 and V01 V0 = function(n,theta) rFJoe(n, 1/theta), V01 = function(V0,theta0,theta1, approx=100000) { alpha <- theta0/theta1 ia <- (V0 > approx) ie <- !ia V01 <- V0 # numeric(length(V0)) V01[ia] <- V0[ia]^(1/alpha) * rstable1(sum(ia),alpha, beta=1, gamma = cos(alpha*pi/2)^(1/alpha), delta = as.integer(alpha==1)) V01[ie] <- sapply(lapply(V0[ie], rFJoe, alpha=alpha),sum) V01 }, ## 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) copJoe@tau(th) - tau, interval = c(1.001, 100), 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) } ) ## {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)