[Vinecopula-commits] r78 - in pkg: . R inst man tests tests/Examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mo Jan 26 16:35:22 CET 2015
Author: etobi
Date: 2015-01-26 16:35:21 +0100 (Mon, 26 Jan 2015)
New Revision: 78
Modified:
pkg/DESCRIPTION
pkg/R/BiCopPar2Tau.r
pkg/R/BiCopTau2Par.r
pkg/inst/ChangeLog
pkg/man/BiCopChiPlot.Rd
pkg/man/BiCopGofTest.Rd
pkg/man/BiCopKPlot.Rd
pkg/man/BiCopLambda.Rd
pkg/man/BiCopPar2Tau.Rd
pkg/man/BiCopSelect.Rd
pkg/man/BiCopTau2Par.Rd
pkg/man/BiCopVuongClarke.Rd
pkg/man/RVineClarkeTest.Rd
pkg/man/RVineGofTest.Rd
pkg/man/RVineMLE.Rd
pkg/man/RVineVuongTest.Rd
pkg/man/VineCopula-package.Rd
pkg/tests/Examples/VineCopula-Ex.Rout.save
pkg/tests/additonalExampleRuns.R
pkg/tests/additonalExampleRuns.Rout.save
Log:
- code cosmetics of additonalExampleRuns.R and the corresponding examples in the manual
- New functionality:
* BiCopTau2Par and BiCopPar2Tau: fully vectorized (parameter/tau input), and sanity checks extendend. Before vector input was not prohibited. However, both functions were not intended to be used for vectorized input.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/DESCRIPTION 2015-01-26 15:35:21 UTC (rev 78)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical Inference of Vine Copulas
Version: 1.4
-Date: 2015-01-20
+Date: 2015-01-26
Author: Ulf Schepsmeier, Jakob Stoeber, Eike Christian Brechmann, Benedikt Graeler
Maintainer: Tobias Erhardt <tobias.erhardt at tum.de>
Depends: R (>= 2.11.0)
Modified: pkg/R/BiCopPar2Tau.r
===================================================================
--- pkg/R/BiCopPar2Tau.r 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/R/BiCopPar2Tau.r 2015-01-26 15:35:21 UTC (rev 78)
@@ -1,209 +1,260 @@
-BiCopPar2Tau<-function(family,par,par2=0)
-{
- if(!(family %in% c(0,1,2,3,4,5,6,7,8,9,10,13,14,16,17,18,19,20,23,24,26,27,28,29,30,33,34,36,37,38,39,40,41,42,51,52,61,62,71,72,104,114,124,134,204,214,224,234))) stop("Copula family not implemented.")
- if(c(7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,42,52,62,72,104,114,124,134,204,214,224,234) %in% family && par2==0) stop("For BB1, BB6, BB7, BB8 and Tawn copulas, 'par2' must be set.")
- if(c(1,3,4,5,6,11,13,14,16,23,24,26,33,34,36,41,51,61,71) %in% family && length(par)<1) stop("'par' not set.")
+BiCopPar2Tau <- function(family, par, par2 = 0){
+
+ ## sanity checks
+ if(length(family) != 1){
+ stop("Input for family has to be a scalar/integer.")
+ }
+ if(!(family %in% c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
+ 13, 14, 16, 17, 18, 19, 20,
+ 23, 24, 26, 27, 28, 29, 30,
+ 33, 34, 36, 37, 38, 39, 40,
+ 41, 42, 51, 52, 61, 62, 71, 72,
+ 104, 114, 124, 134,
+ 204, 214, 224, 234)
+ )){
+ stop("Copula family not implemented.")
+ }
+
+ if(missing(par)){
+ stop("'par' not set.")
+ }
+
+ if(length(par2) > 1 && length(par) != length(par2)){
+ stop("Input for 'par' and 'par2 has to be vectors of same length.")
+ }
+
+ if(family %in% c(7, 8, 9, 10,
+ 17, 18, 19, 20,
+ 27, 28, 29, 30,
+ 37, 38, 39, 40,
+ 42, 52, 62, 72,
+ 104, 114, 124, 134,
+ 204, 214, 224, 234)
+ && par2 == 0){
+ stop("For BB1, BB6, BB7, BB8 and Tawn copulas, 'par2' must be set.")
+ }
- if((family==1 || family==2) && abs(par[1])>=1) stop("The parameter of the Gaussian and t-copula has to be in the interval (-1,1).")
+ if((family==1 || family==2) && any(abs(par[1])>=1)){
+ stop("The (first) parameter of the Gaussian and t-copula has to be in the interval (-1,1).")
+ }
#if(family==2 && par2<=2) stop("The degrees of freedom parameter of the t-copula has to be larger than 2.")
- if((family==3 || family==13) && par<=0) stop("The parameter of the Clayton copula has to be positive.")
- if((family==4 || family==14) && par<1) stop("The parameter of the Gumbel copula has to be in the interval [1,oo).")
- if((family==6 || family==16) && par<=1) stop("The parameter of the Joe copula has to be in the interval (1,oo).")
- if(family==5 && par==0) stop("The parameter of the Frank copula has to be unequal to 0.")
- if((family==7 || family==17) && par<=0) stop("The first parameter of the BB1 copula has to be positive.")
- if((family==7 || family==17) && par2<1) stop("The second parameter of the BB1 copula has to be in the interval [1,oo).")
- if((family==8 || family==18) && par<=0) stop("The first parameter of the BB6 copula has to be in the interval [1,oo).")
- if((family==8 || family==18) && par2<1) stop("The second parameter of the BB6 copula has to be in the interval [1,oo).")
- if((family==9 || family==19) && par<1) stop("The first parameter of the BB7 copula has to be in the interval [1,oo).")
- if((family==9 || family==19) && par2<=0) stop("The second parameter of the BB7 copula has to be positive.")
- if((family==10 || family==20) && par<1) stop("The first parameter of the BB8 copula has to be in the interval [1,oo).")
- if((family==10 || family==20) && (par2<=0 || par2>1)) stop("The second parameter of the BB8 copula has to be in the interval (0,1].")
- if((family==23 || family==33) && par>=0) stop("The parameter of the rotated Clayton copula has to be negative.")
- if((family==24 || family==34) && par>-1) stop("The parameter of the rotated Gumbel copula has to be in the interval (-oo,-1].")
- if((family==26 || family==36) && par>=-1) stop("The parameter of the rotated Joe copula has to be in the interval (-oo,-1).")
- if((family==27 || family==37) && par>=0) stop("The first parameter of the rotated BB1 copula has to be negative.")
- if((family==27 || family==37) && par2>-1) stop("The second parameter of the rotated BB1 copula has to be in the interval (-oo,-1].")
- if((family==28 || family==38) && par>=0) stop("The first parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
- if((family==28 || family==38) && par2>-1) stop("The second parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
- if((family==29 || family==39) && par>-1) stop("The first parameter of the rotated BB7 copula has to be in the interval (-oo,-1].")
- if((family==29 || family==39) && par2>=0) stop("The second parameter of the rotated BB7 copula has to be negative.")
- if((family==30 || family==40) && par>-1) stop("The first parameter of the rotated BB8 copula has to be in the interval (-oo,-1].")
- if((family==30 || family==40) && (par2>=0 || par2<(-1))) stop("The second parameter of the rotated BB8 copula has to be in the interval [-1,0).")
- if((family==41 || family==51) && par<=0) stop("The parameter of the reflection asymmetric copula has to be positive.")
- if((family==61 || family==71) && par>=0) stop("The parameter of the rotated reflection asymmetric copula has to be negative.")
- if(family==42)
- {
- a=par
- b=par2
- limA=(b - 3 - sqrt(9 + 6 * b - 3 * b^2))/2
- if(abs(b)>1) stop("The second parameter of the two-parametric asymmetric copulas has to be in the interval [-1,1]")
- if(a>1 || a<limA) stop("The first parameter of the two-parametric asymmetric copula has to be in the interval [limA(par2),1]")
+ if((family==3 || family==13) && any(par<=0)){
+ stop("The parameter of the (survival) Clayton copula has to be positive.")
}
- if ((family==104 || family==114 || family==204 || family==214) && par<1) stop("Please choose 'par' of the Tawn copula in [1,oo).")
- if ((family==104 || family==114 || family==204 || family==214) && (par2<0 || par2>1)) stop("Please choose 'par2' of the Tawn copula in [0,1].")
- if ((family==124 || family==134 || family==224 || family==234) && par>-1) stop("Please choose 'par' of the Tawn copula in (-oo,-1].")
- if ((family==124 || family==134 || family==224 || family==234) && (par2<0 || par2>1)) stop("Please choose 'par2' of the Tawn copula in [0,1].")
-
- if(family==0)
- {
- tau=0
+ if((family==4 || family==14) && any(par<1)){
+ stop("The parameter of the (survival) Gumbel copula has to be in the interval [1,oo).")
}
- else if(family==1 | family==2)
- {
- tau=2/pi*asin(par)
+ if((family==6 || family==16) && any(par<=1)){
+ stop("The parameter of the (survival) Joe copula has to be in the interval (1,oo).")
}
- else if(family==3 || family==13)
- {
- tau=par/(par+2)
+ if(family==5 && any(par==0)){
+ stop("The parameter of the Frank copula has to be different from 0.")
}
- else if(family==4 || family==14)
- {
- tau=1-1/par
+ if((family==7 || family==17) && any(par<=0)){
+ stop("The first parameter of the BB1 copula has to be positive.")
}
- else if(family==5)
- {
- f=function(x) {x/(exp(x)-1)}
- if(par>0)
- {
- tau=1-4/par+4/par^2*integrate(f,lower=0, upper=par)$value
- }
- else
- {
- tau=1-4/par-4/par^2*integrate(f,lower=par, upper=0)$value
- }
+ if((family==7 || family==17) && any(par2<1)){
+ stop("The second parameter of the BB1 copula has to be in the interval [1,oo).")
}
- else if(family==6 || family==16)
- {
- #tau=1+4/par^2*integrate(function(x) log(x)*x*(1-x)^(2*(1-par)/par), 0, 1)$value
- param1=2/par+1
- tem=digamma(2)-digamma(param1)
- tau=1+tem*2/(2-par)
- tau[par==2]=1-trigamma(2)
+ if((family==8 || family==18) && any(par<=0)){
+ stop("The first parameter of the BB6 copula has to be in the interval [1,oo).")
}
- else if(family==7 || family==17)
- {
- theta=par
- delta=par2
- tau=1-2/(delta*(theta+2))
+ if((family==8 || family==18) && any(par2<1)){
+ stop("The second parameter of the BB6 copula has to be in the interval [1,oo).")
}
- else if(family==8 || family==18)
- {
- theta=par
- delta=par2
- kt<-function(t) {- log(-(1-t)^theta+1)*(1-t-(1-t)^(-theta)+(1-t)^(-theta)*t)/(delta*theta)}
- tau=1+4*integrate(kt,0,1)$value
+ if((family==9 || family==19) && any(par<1)){
+ stop("The first parameter of the BB7 copula has to be in the interval [1,oo).")
}
- else if(family==9 || family==19)
- {
- #theta=par
- #delta=par2
- #tau=1-2/(delta*(2-theta))+4/(theta^2*delta)*gamma(delta+2)*gamma((2-2*theta)/(theta)+1)/gamma(delta+3+(2-2*theta)/(theta))
- kt<-function(t) {( (1-(1-t)^par)^-par2-1 )/( -par*par2*(1-t)^(par-1)*(1-(1-t)^par)^(-par2-1) )}
- tau=1+4*integrate(kt,0,1)$value
- #kt<-function(t) { 1/( (1-t)^(par-1) ) }
- #kt2<-function(t) { 1-t }
- #kt3<-function(t) { 1/( (1-t)^(par-1)*(1-(1-t)^par)^(-par2-1) ) }
- #tau=1-4/par/par2*(integrate(kt,0,1)$value-integrate(kt2,0,1)$value-integrate(kt3,0,1)$value)
-
+ if((family==9 || family==19) && any(par2<=0)){
+ stop("The second parameter of the BB7 copula has to be positive.")
}
- else if(family==10 || family==20)
- {
- theta=par
- delta=par2
- kt<-function(t) {-log(((1-t*delta)^theta-1)/((1-delta)^theta-1))*(1-t*delta-(1-t*delta)^(-theta)+(1-t*delta)^(-theta)*t*delta)/(theta*delta)}
- tau=1+4*integrate(kt,0,1)$value
+ if((family==10 || family==20) && any(par<1)){
+ stop("The first parameter of the BB8 copula has to be in the interval [1,oo).")
}
- else if(family==23 || family==33)
- {
- tau=par/(-par+2)
+ if((family==10 || family==20) && any(par2<=0 || par2>1)){
+ stop("The second parameter of the BB8 copula has to be in the interval (0,1].")
}
- else if(family==24 || family==34)
- {
- tau=-1-1/par
+ if((family==23 || family==33) && any(par>=0)){
+ stop("The parameter of the (90/270 degree) rotated Clayton copula has to be negative.")
}
- else if(family==26 || family==36)
- {
- #tau=-1-4/par^2*integrate(function(x) log(x)*x*(1-x)^(2*(1+par)/-par), 0, 1)$value
- theta=-par
- param1=2/theta+1
- tem=digamma(2)-digamma(param1)
- tau=1+tem*2/(2-theta)
- tau[theta==2]=1-trigamma(2)
- tau=-tau
+ if((family==24 || family==34) && any(par>-1)){
+ stop("The parameter of the (90/270 degree) rotated Gumbel copula has to be in the interval (-oo,-1].")
}
- else if(family==27 || family==37)
- {
- theta=-par
- delta=-par2
- tau=1-2/(delta*(theta+2))
- tau=-tau
+ if((family==26 || family==36) && any(par>=-1)){
+ stop("The parameter of the (90/270 degree) rotated Joe copula has to be in the interval (-oo,-1).")
}
- else if(family==28 || family==38)
- {
- theta=-par
- delta=-par2
- kt<-function(t) {- log(-(1-t)^theta+1)*(1-t-(1-t)^(-theta)+(1-t)^(-theta)*t)/(delta*theta)}
- tau=1+4*integrate(kt,0,1)$value
- tau=-tau
+ if((family==27 || family==37) && any(par>=0)){
+ stop("The first parameter of the (90/270 degree) rotated BB1 copula has to be negative.")
}
- else if(family==29 || family==39)
- {
- theta=-par
- delta=-par2
- #tau=1-2/(delta*(2-theta))+4/(theta^2*delta)*gamma(delta+2)*gamma((2-2*theta)/(theta)+1)/gamma(delta+3+(2-2*theta)/(theta))
- kt<-function(t) {( (1-(1-t)^theta)^(-delta)-1 )/( -theta*delta*(1-t)^(theta-1)*(1-(1-t)^theta)^(-delta-1) )}
- tau=1+4*integrate(kt,0,1)$value
- tau=-tau
+ if((family==27 || family==37) && any(par2>-1)){
+ stop("The second parameter of the (90/270 degree) rotated BB1 copula has to be in the interval (-oo,-1].")
}
- else if(family==30 || family==40)
- {
- theta=-par
- delta=-par2
- kt<-function(t) {-log(((1-t*delta)^theta-1)/((1-delta)^theta-1))*(1-t*delta-(1-t*delta)^(-theta)+(1-t*delta)^(-theta)*t*delta)/(theta*delta)}
- tau=1+4*integrate(kt,0,1)$value
- tau=-tau
+ if((family==28 || family==38) && any(par>=0)){
+ stop("The first parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
}
- else if(family==41 || family==51)
- {
- de=par
- ln2=log(2)
- tem=(2-2*de)*ln2+lgamma(2*de)-2*lgamma(1+de)
- tau=1-de*exp(tem)
+ if((family==28 || family==38) && any(par2>-1)){
+ stop("The second parameter of the (90/270 degree) rotated BB6 copula has to be in the interval (-oo,-1].")
}
- else if(family==61 || family==71)
- {
- de=-par
- ln2=log(2)
- tem=(2-2*de)*ln2+lgamma(2*de)-2*lgamma(1+de)
- tau=1-de*exp(tem)
- tau=-tau
+ if((family==29 || family==39) && any(par>-1)){
+ stop("The first parameter of the (90/270 degree) rotated BB7 copula has to be in the interval (-oo,-1].")
}
- else if(family==42)
- {
- tau=(75*par2-par2^2+par*(25-par2))/450
+ if((family==29 || family==39) && any(par2>=0)){
+ stop("The second parameter of the (90/270 degree) rotated BB7 copula has to be negative.")
}
- else if(family==104 || family==114 || family==204 || family==214)
- {
- par3=1
- tau_int=function(t)
- {
- Afunc = .C("Tawn2",as.double(t),as.integer(length(t)),as.double(par),as.double(par2),as.double(1),as.double(rep(0,length(t))),PACKAGE='VineCopula')[[6]]
- Afunc2Deriv = .C("d2Tawn",as.double(t),as.integer(length(t)),as.double(par),as.double(par2),as.double(1),as.double(rep(0,length(t))),PACKAGE='VineCopula')[[6]]
- (t*(1-t))*Afunc2Deriv/ Afunc
- }
- tau<-integrate(tau_int,0,1)[[1]]
+ if((family==30 || family==40) && any(par>-1)){
+ stop("The first parameter of the (90/270 degree) rotated BB8 copula has to be in the interval (-oo,-1].")
}
- else if(family==124 || family==134 || family==224 || family==234)
- {
- par3=1
- tau_int=function(t)
- {
- Afunc = .C("Tawn2",as.double(t),as.integer(length(t)),as.double(-par),as.double(par2),as.double(1),as.double(rep(0,length(t))),PACKAGE='VineCopula')[[6]]
- Afunc2Deriv = .C("d2Tawn",as.double(t),as.integer(length(t)),as.double(-par),as.double(par2),as.double(1),as.double(rep(0,length(t))),PACKAGE='VineCopula')[[6]]
- (t*(1-t))*Afunc2Deriv/ Afunc
+ if((family==30 || family==40) && any(par2>=0 || par2<(-1))){
+ stop("The second parameter of the (90/270 degree) rotated BB8 copula has to be in the interval [-1,0).")
+ }
+ if((family==41 || family==51) && any(par<=0)){
+ stop("The parameter of the (survival) reflection asymmetric copula has to be positive.")
+ }
+ if((family==61 || family==71) && any(par>=0)){
+ stop("The parameter of the (90/270 degree) rotated reflection asymmetric copula has to be negative.")
+ }
+ if(family==42){
+ a <- par
+ b <- par2
+ limA <- (b - 3 - sqrt(9 + 6 * b - 3 * b^2))/2
+ if(any(abs(b)>1)) stop("The second parameter of the two-parametric asymmetric copulas has to be in the interval [-1,1]")
+ if(any(a>1 || a<limA)) stop("The first parameter of the two-parametric asymmetric copula has to be in the interval [limA(par2),1]")
+ }
+ if((family==104 || family==114 || family==204 || family==214) && any(par<1)){
+ stop("Please choose 'par' of the (180 degree rotated) Tawn copula in [1,oo).")
+ }
+ if((family==104 || family==114 || family==204 || family==214) && any(par2<0 || par2>1)){
+ stop("Please choose 'par2' of the (180 degree rotated) Tawn copula in [0,1].")
+ }
+ if((family==124 || family==134 || family==224 || family==234) && any(par>-1)){
+ stop("Please choose 'par' of the (90/270 degree) rotated Tawn copula in (-oo,-1].")
+ }
+ if((family==124 || family==134 || family==224 || family==234) && any(par2<0 || par2>1)){
+ stop("Please choose 'par2' of the (90/270 degree) rotated Tawn copula in [0,1].")
+ }
+
+ ## calculation of tau(s) depending on pair-copula family
+ if(family==0){
+ tau <- rep(0, times = length(par))
+ } else if(family==1 | family==2){
+ tau <- 2/pi*asin(par)
+ } else if(family==3 || family==13){
+ tau <- par/(par+2)
+ } else if(family==4 || family==14){
+ tau <- 1-1/par
+ } else if(family==5){
+ f <- function(x) x/(exp(x)-1)
+ fu <- function(x) integrate(f, lower = 0, upper = x)$value
+ fl <- function(x) integrate(f, lower = x, upper = 0)$value
+ if(any(par>0)){
+ tau <- 1 - 4/par + 4/par^2 * sapply(par, fu)
+ } else {
+ tau <- 1 - 4/par - 4/par^2 * sapply(par, fl)
}
- tau<-integrate(tau_int,0,1)[[1]]
- tau=-tau
+ } else if(family==6 || family==16){
+ # tau = 1 + 4/par^2 * integrate(function(x) log(x)*x*(1-x)^(2*(1-par)/par), 0, 1)$value
+ param1 <- 2/par + 1
+ tem <- digamma(2) - digamma(param1)
+ tau <- 1 + tem*2/(2 - par)
+ tau[par==2] <- 1 - trigamma(2)
+ } else if(family==7 || family==17){
+ theta <- par
+ delta <- par2
+ tau <- 1 - 2/(delta*(theta+2))
+ } else if(family==8 || family==18){
+ theta <- par
+ delta <- par2
+ kt <- function(t, th, de){-log(-(1-t)^th+1)*(1-t-(1-t)^(-th)+(1-t)^(-th)*t)/(de*th)}
+ tau <- 1 + 4*mapply(function(theta, delta){integrate(function(t){kt(t, th=theta, de=delta)}, 0, 1)$value},
+ theta, delta)
+ } else if(family==9 || family==19){
+ theta <- par
+ delta <- par2
+ #tau=1-2/(delta*(2-theta))+4/(theta^2*delta)*gamma(delta+2)*gamma((2-2*theta)/(theta)+1)/gamma(delta+3+(2-2*theta)/(theta))
+ kt <- function(t, th, de){( (1-(1-t)^th)^-de-1 )/( -th*de*(1-t)^(th-1)*(1-(1-t)^th)^(-de-1) )}
+ tau <- 1 + 4*mapply(function(theta, delta){integrate(function(t){kt(t, th=theta, de=delta)}, 0, 1)$value},
+ theta, delta)
+ #kt <- function(t) { 1/( (1-t)^(par-1) ) }
+ #kt2 <- function(t) { 1-t }
+ #kt3 <- function(t) { 1/( (1-t)^(par-1)*(1-(1-t)^par)^(-par2-1) ) }
+ #tau <- 1-4/par/par2*(integrate(kt,0,1)$value-integrate(kt2,0,1)$value-integrate(kt3,0,1)$value)
+ } else if(family==10 || family==20){
+ theta <- par
+ delta <- par2
+ kt <- function(t, th, de){-log(((1-t*de)^th-1)/((1-de)^th-1))*(1-t*de-(1-t*de)^(-th)+(1-t*de)^(-th)*t*de)/(th*de)}
+ tau <- 1 + 4*mapply(function(theta, delta){integrate(function(t){kt(t, th=theta, de=delta)}, 0, 1)$value},
+ theta, delta)
+ } else if(family==23 || family==33){
+ tau <- par/(-par+2)
+ } else if(family==24 || family==34){
+ tau <- -1-1/par
+ } else if(family==26 || family==36){
+ #tau <- -1-4/par^2*integrate(function(x) log(x)*x*(1-x)^(2*(1+par)/-par), 0, 1)$value
+ theta <- -par
+ param1 <- 2/theta+1
+ tem <- digamma(2)-digamma(param1)
+ tau <- 1+tem*2/(2-theta)
+ tau[theta==2] <- 1-trigamma(2)
+ tau <- -tau
+ } else if(family==27 || family==37){
+ theta <- -par
+ delta <- -par2
+ tau <- 1-2/(delta*(theta+2))
+ tau <- -tau
+ } else if(family==28 || family==38){
+ theta <- -par
+ delta <- -par2
+ kt <- function(t, th, de){- log(-(1-t)^th+1)*(1-t-(1-t)^(-th)+(1-t)^(-th)*t)/(de*th)}
+ tau <- 1 + 4*mapply(function(theta, delta){integrate(function(t){kt(t, th=theta, de=delta)}, 0, 1)$value},
+ theta, delta)
+ tau <- -tau
+ } else if(family==29 || family==39){
+ theta <- -par
+ delta <- -par2
+ #tau <- 1-2/(delta*(2-theta))+4/(theta^2*delta)*gamma(delta+2)*gamma((2-2*theta)/(theta)+1)/gamma(delta+3+(2-2*theta)/(theta))
+ kt <- function(t, th, de){( (1-(1-t)^th)^(-de)-1 )/( -th*de*(1-t)^(th-1)*(1-(1-t)^th)^(-de-1) )}
+ tau <- 1 + 4*mapply(function(theta, delta){integrate(function(t){kt(t, th=theta, de=delta)}, 0, 1)$value},
+ theta, delta)
+ tau <- -tau
+ } else if(family==30 || family==40){
+ theta <- -par
+ delta <- -par2
+ kt <- function(t, th, de){-log(((1-t*de)^th-1)/((1-de)^th-1))*(1-t*de-(1-t*de)^(-th)+(1-t*de)^(-th)*t*de)/(th*de)}
+ tau <- 1 + 4*mapply(function(theta, delta){integrate(function(t){kt(t, th=theta, de=delta)}, 0, 1)$value},
+ theta, delta)
+ tau <- -tau
+ } else if(family==41 || family==51){
+ de <- par
+ ln2 <- log(2)
+ tem <- (2 - 2*de)*ln2 + lgamma(2*de) - 2*lgamma(1 + de)
+ tau <- 1 - de*exp(tem)
+ } else if(family==61 || family==71){
+ de <- -par
+ ln2 <- log(2)
+ tem <- (2 - 2*de)*ln2 + lgamma(2*de) - 2*lgamma(1 + de)
+ tau <- 1 - de*exp(tem)
+ tau <- -tau
+ } else if(family==42){
+ tau <- (75*par2 - par2^2 + par*(25 - par2))/450
+ } else if(family==104 || family==114 || family==204 || family==214){
+ par3 <- 1
+ tau_int <- function(t, th, de){
+ Afunc = .C("Tawn2", as.double(t), as.integer(length(t)), as.double(th), as.double(de), as.double(1), as.double(rep(0,length(t))), PACKAGE='VineCopula')[[6]]
+ Afunc2Deriv = .C("d2Tawn", as.double(t), as.integer(length(t)), as.double(th), as.double(de), as.double(1), as.double(rep(0,length(t))), PACKAGE='VineCopula')[[6]]
+ (t*(1 - t)) * Afunc2Deriv/Afunc
+ }
+ tau <- mapply(function(par, par2){integrate(function(t){tau_int(t, th=par, de=par2)}, 0, 1)$value},
+ par, par2)
+ } else if(family==124 || family==134 || family==224 || family==234){
+ par3 <- 1
+ tau_int <- function(t, th, de){
+ Afunc = .C("Tawn2", as.double(t), as.integer(length(t)), as.double(-th), as.double(de), as.double(1), as.double(rep(0,length(t))), PACKAGE='VineCopula')[[6]]
+ Afunc2Deriv = .C("d2Tawn", as.double(t), as.integer(length(t)), as.double(-th), as.double(de), as.double(1), as.double(rep(0,length(t))), PACKAGE='VineCopula')[[6]]
+ (t*(1 - t)) * Afunc2Deriv/Afunc
+ }
+ tau <- mapply(function(par, par2){integrate(function(t){tau_int(t, th=par, de=par2)}, 0, 1)$value},
+ par, par2)
+ tau <- -tau
}
return(tau)
Modified: pkg/R/BiCopTau2Par.r
===================================================================
--- pkg/R/BiCopTau2Par.r 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/R/BiCopTau2Par.r 2015-01-26 15:35:21 UTC (rev 78)
@@ -1,38 +1,50 @@
-BiCopTau2Par<-function(family,tau){
+BiCopTau2Par <- function(family, tau){
+
+ ## sanity checks
+ if(length(family) != 1){
+ stop("Input for family has to be a scalar/integer.")
+ }
+ if(!(family %in% c(0, 1, 2, 3, 4, 5, 6,
+ 13, 14, 16,
+ 23, 24, 26,
+ 33, 34, 36,
+ 41, 51, 61, 71)
+ )){
+ stop("Copula family not implemented.")
+ }
- if(!(family %in% c(0,1,3,4,5,6,13,14,16,23,24,26,33,34,36,41,51,61,71))) stop("Copula family not implemented.")
-
+ ## calculation of parameter(s) depending on pair-copula family
if(family == 0){
- par = 0
- }else if(family == 1){
- par = sin(pi*tau/2)
+ par <- rep(0, times = length(tau))
+ }else if(family %in% 1:2){
+ par <- sin(pi*tau/2)
}else if(family %in% c(3,13)){
- if(tau<=0) stop("Clayton copula cannot be used for tau<=0.")
- par = 2*tau/(1-tau)
+ if(any(tau<=0)) stop("Clayton copula cannot be used for tau<=0.")
+ par <- 2*tau/(1-tau)
}else if(family %in% c(4,14)){
- if(tau<0) stop("Gumbel copula cannot be used for tau<0.")
- par = 1/(1-tau)
+ if(any(tau<0)) stop("Gumbel copula cannot be used for tau<0.")
+ par <- 1/(1-tau)
}else if(family == 5){
- if(tau==0) stop("Frank copula cannot be used for tau=0.")
- par = Frank.itau.JJ(tau)
+ if(any(tau==0)) stop("Frank copula cannot be used for tau=0.")
+ par <- sapply(tau, Frank.itau.JJ)
}else if(family %in% c(6,16)){
- if(tau<=0) stop("Joe copula cannot be used for tau<=0.")
- par = Joe.itau.JJ(tau)
+ if(any(tau<=0)) stop("Joe copula cannot be used for tau<=0.")
+ par <- sapply(tau, Joe.itau.JJ)
}else if(family %in% c(23,33)){
- if(tau>=0) stop("Rotated Clayton copula cannot be used for tau>=0.")
- par = 2*tau/(1+tau)
+ if(any(tau>=0)) stop("Rotated Clayton copula cannot be used for tau>=0.")
+ par <- 2*tau/(1+tau)
}else if(family %in% c(24,34)){
- if(tau>0) stop("Rotated Gumbel copula cannot be used for tau>0.")
- par = -(1/(1+tau))
+ if(any(tau>0)) stop("Rotated Gumbel copula cannot be used for tau>0.")
+ par <- -(1/(1+tau))
}else if(family %in% c(26,36)){
- if(tau>=0) stop("Rotated Joe copula cannot be used for tau>=0.")
- par = -Joe.itau.JJ(-tau)
+ if(any(tau>=0)) stop("Rotated Joe copula cannot be used for tau>=0.")
+ par <- -sapply(-tau, Joe.itau.JJ)
}
else if(family %in% c(41,51)){
- par = ipsA.tau2cpar(tau)
+ par <- sapply(tau, ipsA.tau2cpar)
}
else if(family %in% c(61,71)){
- par=-ipsA.tau2cpar(-tau)
+ par <- -sapply(-tau, ipsA.tau2cpar)
}
return(par)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/inst/ChangeLog 2015-01-26 15:35:21 UTC (rev 78)
@@ -1,10 +1,15 @@
Changes for R-package VineCopula
-Current authors: Ulf Schepsmeier, Thomas Nagler and Benedikt Graeler
+Current authors: Ulf Schepsmeier, Thomas Nagler, Tobias Erhardt and Benedikt Graeler
Former authors: Eike Brechmann and Jakob Stoeber
Maintainer: Tobias Erhardt <tobias.erhardt at tum.de> and Thomas Nagler <thomas.nagler at tum.de>
+Version 1.4 (January 26, 2015)
+
+- New functionality:
+ * BiCopTau2Par and BiCopPar2Tau: fully vectorized (parameter/tau input), and sanity checks extendend. Before vector input was not prohibited. However, both functions were not intended to be used for vectorized input.
+
Version 1.3-2 (January 19, 2015)
- New author: Thomas Nagler
Modified: pkg/man/BiCopChiPlot.Rd
===================================================================
--- pkg/man/BiCopChiPlot.Rd 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/man/BiCopChiPlot.Rd 2015-01-26 15:35:21 UTC (rev 78)
@@ -77,22 +77,22 @@
\examples{
\dontrun{
# chi-plots for bivariate Gaussian copula data
-n = 500
-tau = 0.5
+n <- 500
+tau <- 0.5
# simulate copula data
-fam = 1
-theta = BiCopTau2Par(fam,tau)
-dat = BiCopSim(n,fam,theta)
+fam <- 1
+theta <- BiCopTau2Par(fam, tau)
+set.seed(123)
+dat <- BiCopSim(n, fam, theta)
# create chi-plots
-dev.new(width=16,height=5)
-par(mfrow=c(1,3))
-BiCopChiPlot(dat[,1],dat[,2],xlim=c(-1,1),ylim=c(-1,1),
+par(mfrow = c(1,3))
+BiCopChiPlot(dat[,1], dat[,2], xlim = c(-1,1), ylim = c(-1,1),
main="General chi-plot")
-BiCopChiPlot(dat[,1],dat[,2],mode="lower",xlim=c(-1,1),
- ylim=c(-1,1),main="Lower chi-plot")
-BiCopChiPlot(dat[,1],dat[,2],mode="upper",xlim=c(-1,1),
- ylim=c(-1,1),main="Upper chi-plot")
+BiCopChiPlot(dat[,1], dat[,2], mode = "lower", xlim = c(-1,1),
+ ylim = c(-1,1), main = "Lower chi-plot")
+BiCopChiPlot(dat[,1], dat[,2], mode = "upper", xlim = c(-1,1),
+ ylim = c(-1,1), main = "Upper chi-plot")
}
}
Modified: pkg/man/BiCopGofTest.Rd
===================================================================
--- pkg/man/BiCopGofTest.Rd 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/man/BiCopGofTest.Rd 2015-01-26 15:35:21 UTC (rev 78)
@@ -124,25 +124,26 @@
\seealso{\code{\link{BiCopDeriv2}}, \code{\link{BiCopDeriv}}, \code{\link{BiCopIndTest}}, \code{\link{BiCopVuongClarke}}}
\examples{
-# simulate from a bivariate Clayton copula
-simdata = BiCopSim(300,3,2)
-u1 = simdata[,1]
-u2 = simdata[,2]
+ # simulate from a bivariate Clayton copula
+ set.seed(123)
+ simdata <- BiCopSim(300, 3, 2)
+ u1 <- simdata[,1]
+ u2 <- simdata[,2]
+
+ # perform White's goodness-of-fit test for the true copula
+ BiCopGofTest(u1, u2, family = 3)
+
+ # perform Kendall's goodness-of-fit test for the Frank copula
+ BiCopGofTest(u1, u2, family = 5)
-# perform White's goodness-of-fit test for the true copula
-BiCopGofTest(u1,u2,family=3)
-
-# perform Kendall's goodness-of-fit test for the Frank copula
-BiCopGofTest(u1,u2,family=5)
-
\dontrun{
# perform Kendall's goodness-of-fit test for the true copula
-gof = BiCopGofTest(u1,u2,family=3,method="kendall")
+gof <- BiCopGofTest(u1, u2, family = 3, method = "kendall")
gof$p.value.CvM
gof$p.value.KS
# perform Kendall's goodness-of-fit test for the Frank copula
-gof = BiCopGofTest(u1,u2,family=5,method="kendall")
+gof <- BiCopGofTest(u1, u2, family = 5, method = "kendall")
gof$p.value.CvM
gof$p.value.KS
}
Modified: pkg/man/BiCopKPlot.Rd
===================================================================
--- pkg/man/BiCopKPlot.Rd 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/man/BiCopKPlot.Rd 2015-01-26 15:35:21 UTC (rev 78)
@@ -62,23 +62,24 @@
\examples{
\dontrun{
# Gaussian and Clayton copulas
-n = 500
-tau = 0.5
+n <- 500
+tau <- 0.5
# simulate from Gaussian copula
-fam1 = 1
-theta1 = BiCopTau2Par(fam1,tau)
-dat1 = BiCopSim(n,fam1,theta1)
+fam1 <- 1
+theta1 <- BiCopTau2Par(fam1, tau)
+set.seed(123)
+dat1 <- BiCopSim(n, fam1, theta1)
# simulate from Clayton copula
-fam2 = 3
-theta2 = BiCopTau2Par(fam2,tau)
-dat2 = BiCopSim(n,fam2,theta2)
+fam2 <- 3
+theta2 <- BiCopTau2Par(fam2, tau)
+set.seed(123)
+dat2 <- BiCopSim(n, fam2, theta2)
# create K-plots
-dev.new(width=10,height=5)
par(mfrow=c(1,2))
-BiCopKPlot(dat1[,1],dat1[,2],main="Gaussian copula")
-BiCopKPlot(dat2[,1],dat2[,2],main="Clayton copula")
+BiCopKPlot(dat1[,1], dat1[,2], main = "Gaussian copula")
+BiCopKPlot(dat2[,1], dat2[,2], main = "Clayton copula")
}
}
Modified: pkg/man/BiCopLambda.Rd
===================================================================
--- pkg/man/BiCopLambda.Rd 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/man/BiCopLambda.Rd 2015-01-26 15:35:21 UTC (rev 78)
@@ -89,34 +89,34 @@
\examples{
\dontrun{
# Clayton and rotated Clayton copulas
-n = 1000
-tau = 0.5
+n <- 1000
+tau <- 0.5
# simulate from Clayton copula
-fam = 3
-theta = BiCopTau2Par(fam,tau)
-dat = BiCopSim(n,fam,theta)
+fam <- 3
+theta <- BiCopTau2Par(fam, tau)
+set.seed(123)
+dat <- BiCopSim(n, fam, theta)
# create lambda-function plots
-dev.new(width=16,height=5)
-par(mfrow=c(1,3))
-BiCopLambda(dat[,1],dat[,2]) # empirical lambda-function
-BiCopLambda(family=fam,par=theta) # theoretical lambda-function
-BiCopLambda(dat[,1],dat[,2],family=fam,par=theta) # both
+par(mfrow = c(1,3))
+BiCopLambda(dat[,1], dat[,2]) # empirical lambda-function
+BiCopLambda(family = fam, par = theta) # theoretical lambda-function
+BiCopLambda(dat[,1], dat[,2], family = fam, par = theta) # both
# simulate from rotated Clayton copula (90 degrees)
-fam = 23
-theta = BiCopTau2Par(fam,-tau)
-dat = BiCopSim(n,fam,theta)
-
+fam <- 23
+theta <- BiCopTau2Par(fam, -tau)
+set.seed(123)
+dat <- BiCopSim(n, fam, theta)
+
# rotate the data to standard Clayton copula data
-rot_dat = 1-dat[,1]
+rot_dat <- 1 - dat[,1]
-dev.new(width=16,height=5)
-par(mfrow=c(1,3))
-BiCopLambda(rot_dat,dat[,2]) # empirical lambda-function
-BiCopLambda(family=3,par=-theta) # theoretical lambda-function
-BiCopLambda(rot_dat,dat[,2],family=3,par=-theta) # both
+par(mfrow = c(1,3))
+BiCopLambda(rot_dat, dat[,2]) # empirical lambda-function
+BiCopLambda(family=3, par = -theta) # theoretical lambda-function
+BiCopLambda(rot_dat, dat[,2], family = 3, par = -theta) # both
}
}
Modified: pkg/man/BiCopPar2Tau.Rd
===================================================================
--- pkg/man/BiCopPar2Tau.Rd 2015-01-20 18:18:31 UTC (rev 77)
+++ pkg/man/BiCopPar2Tau.Rd 2015-01-26 15:35:21 UTC (rev 78)
@@ -54,18 +54,18 @@
\code{224} = rotated Tawn type 2 copula (90 degrees) \cr
\code{234} = rotated Tawn type 2 copula (270 degrees) \cr
}
- \item{par}{Copula parameter.}
- \item{par2}{Second parameter for the two parameter t-, BB1, BB6, BB7, BB8, Tawn type 1 and type 2 copulas (default: \code{par2 = 0}).
+ \item{par}{Copula parameter (vector).}
+ \item{par2}{Second parameter (vector of same length as \code{par}) for the two parameter t-, BB1, BB6, BB7, BB8, Tawn type 1 and type 2 copulas (default: \code{par2 = 0}).
Note that the degrees of freedom parameter of the t-copula does not need to be set,
because the theoretical Kendall's tau value of the t-copula is independent of this choice.}
}
\value{
-Theoretical value of Kendall's tau corresponding to the bivariate copula family and parameter(s)
+Theoretical value of Kendall's tau (vector) corresponding to the bivariate copula family and parameter(vectors)
(\eqn{\theta} for one parameter families and the first parameter of the t-copula,
-\eqn{\theta} and \eqn{\delta} for the two parameter BB1, BB6, BB7 and BB8 copulas).
+\eqn{\theta} and \eqn{\delta} for the two parameter BB1, BB6, BB7, BB8, Tawn type 1 and type 2 copulas).
\tabular{ll}{
-No. \tab Kendall's tau \cr
+No. (\code{family}) \tab Kendall's tau (\code{tau}) \cr
\code{1, 2} \tab \eqn{\frac{2}{\pi}\arcsin(\theta)}{2 / \pi arcsin(\theta)} \cr
\code{3, 13} \tab \eqn{\frac{\theta}{\theta+2}}{\theta / (\theta+2)} \cr
\code{4, 14} \tab \eqn{1-\frac{1}{\theta}}{1-1/\theta} \cr
@@ -97,7 +97,7 @@
}
}
-\author{Ulf Schepsmeier}
+\author{Ulf Schepsmeier, Tobias Erhardt}
\references{
Joe, H. (1997).
@@ -113,13 +113,28 @@
\examples{
## Example 1: Gaussian copula
-tt1 = BiCopPar2Tau(1,0.7)
+tau0 <- 0.5
+rho <- BiCopTau2Par(family = 1, tau = tau0)
# transform back
-BiCopTau2Par(1,tt1)
+tau <- BiCopPar2Tau(family = 1, par = rho)
+tau - 2/pi*asin(rho)
## Example 2: Clayton copula
-BiCopPar2Tau(3,1.3)
+theta <- BiCopTau2Par(family = 3, tau = c(0.4,0.5,0.6))
+BiCopPar2Tau(family = 3, par = theta)
+
+
+## Example 3:
+vpar <- seq(from = 1.1, to = 10, length.out = 100)
+tauC <- BiCopPar2Tau(family = 3, par = vpar)
+tauG <- BiCopPar2Tau(family = 4, par = vpar)
+tauF <- BiCopPar2Tau(family = 5, par = vpar)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 78
Mehr Informationen über die Mailingliste Vinecopula-commits