[Vinecopula-commits] r25 - in pkg: . R inst man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mi Okt 9 15:49:17 CEST 2013
Author: ulf
Date: 2013-10-09 15:49:16 +0200 (Wed, 09 Oct 2013)
New Revision: 25
Added:
pkg/R/AD.R
pkg/R/BetaMatrix.r
pkg/R/BiCopPar2Beta.r
pkg/R/ChatZj.R
pkg/R/CvM.R
pkg/R/Fhat.R
pkg/R/KS.R
pkg/R/RVineGofTest3.r
pkg/R/RVinePIT.r
pkg/R/RVinePar2Beta.r
pkg/R/gof_ECP.r
pkg/R/gof_PIT.r
pkg/R/gof_White.r
pkg/man/BetaMatrix.Rd
pkg/man/BiCopPar2Beta.Rd
pkg/man/RVineGofTest.Rd
pkg/man/RVinePIT.Rd
pkg/man/RVinePar2Beta.Rd
Removed:
pkg/VineCopula.pdf
Modified:
pkg/R/BiCopCDF.r
pkg/R/BiCopEst.r
pkg/R/BiCopGofTest.r
pkg/R/BiCopHfunc.r
pkg/R/BiCopMetaContour.r
pkg/R/BiCopName.r
pkg/R/BiCopPDF.r
pkg/R/BiCopPar2TailDep.r
pkg/R/BiCopPar2Tau.r
pkg/R/BiCopSelect.r
pkg/R/BiCopSim.R
pkg/R/BiCopVuongClarke.r
pkg/R/RVineCopSelect.r
pkg/R/RVineMLE.R
pkg/R/RVineMatrix.R
pkg/inst/ChangeLog
pkg/man/BiCopGofTest.Rd
pkg/man/VineCopula-package.Rd
pkg/src/hfunc.c
pkg/src/likelihood.c
Log:
Grosses update. Neuerungen s. ChangeLog
Was jetzt noch fehlt ist die Dokumentation der Tawn copula. Die muss ich noch in alle Hilfefiles reinschreiben.
Added: pkg/R/AD.R
===================================================================
--- pkg/R/AD.R (rev 0)
+++ pkg/R/AD.R 2013-10-09 13:49:16 UTC (rev 25)
@@ -0,0 +1,20 @@
+"AD" =
+function(cdf=NULL)
+{
+ # Cumulative distribution function test:
+ # Function that computes the Anderson-Darling test statistic
+ #--------------------------------------------------------------------------
+ # INPUT:
+ # cdf CDF for which to compute AD test
+ # OUTPUT:
+ # AD Anderson-Darling test statistic
+ #--------------------------------------------------------------------------
+ # Author: Daniel Berg <daniel at danielberg.no>
+ # Date: 27.03.2006
+ # Version: 1.0.1
+ #--------------------------------------------------------------------------
+ n = length(cdf)
+ AD = .C("ADtest",as.double(cdf),as.integer(n),as.double(0),PACKAGE='VineCopula')[[3]]
+ AD
+}
+
Added: pkg/R/BetaMatrix.r
===================================================================
--- pkg/R/BetaMatrix.r (rev 0)
+++ pkg/R/BetaMatrix.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -0,0 +1,48 @@
+BetaMatrix<-function(data)
+{
+ d<-dim(data)[2]
+
+ betahat=matrix(1,d,d)
+ for(i in 1:(d-1))
+ {
+ u1=data[,i]
+ for(j in (i+1):d)
+ {
+ u2=data[,j]
+ betahat[i,j]<-betaFunc(u1,u2,1/2,1/2)
+ betahat[j,i]=betahat[i,j]
+ }
+ }
+
+return(betahat)
+}
+
+
+# empirical copula
+empcop<-function(u1,u2,u,v)
+{
+ n=length(u1)
+ a<-which(u1<u)
+ b<-which(u2<v)
+ sc<-intersect(a,b)
+ return(1/n*length(sc))
+}
+
+# survival copula
+survivalcop<-function(u1,u2,u,v)
+{
+ n=length(u1)
+ a<-which(u1>u)
+ b<-which(u2>v)
+ sc<-intersect(a,b)
+ return(1/n*length(sc))
+}
+
+# h_d
+h<-function(u,v) (min(u,v)+min(1-u)-u*v-(1-u)*(1-v))^-1
+
+# g_d
+g<-function(u,v) (u*v)+(1-u)*(1-v)
+
+# beta
+betaFunc<-function(u1,u2,u,v) h(u,v)*(empcop(u1,u2,u,v)+survivalcop(u1,u2,u,v)-g(u,v))
\ No newline at end of file
Modified: pkg/R/BiCopCDF.r
===================================================================
--- pkg/R/BiCopCDF.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopCDF.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -5,8 +5,8 @@
if(any(u2>1) || any(u2<0)) stop("Data has be in the interval [0,1].")
if(length(u1)!=length(u2)) stop("Lengths of 'u1' and 'u2' do not match.")
if(family==2) stop("The CDF of the t-copula is not implemented.")
- if(!(family %in% c(0,1,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,51,61,71))) stop("Copula family not implemented.")
- if(family %in% c(7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40) && par2==0) stop("For BB1, BB6, BB7 and BB8 copulas, 'par2' must be set.")
+ if(!(family %in% c(0,1,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,51,61,71,104,114,124,134,204,214,224,234))) stop("Copula family not implemented.")
+ if(family %in% c(7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,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 %in% c(1,3,4,5,6,13,14,16,23,24,26,33,34,36,41,51,61,71) && length(par)<1) stop("'par' not set.")
if((family==1) && abs(par[1])>=1) stop("The parameter of the Gaussian has to be in the interval (-1,1).")
@@ -36,6 +36,10 @@
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==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].")
res = rep(NA, length(u1))
@@ -56,6 +60,59 @@
res = u2-.C("archCDF",as.double(1-u1),as.double(u2),as.integer(length(u1)),as.double(c(-par,-par2)),as.integer(family-20),as.double(rep(0,length(u1))),PACKAGE='VineCopula')[[6]]
}else if(family %in% c(33,34,36:40,71)){
res = u1-.C("archCDF",as.double(u1),as.double(1-u2),as.integer(length(u1)),as.double(c(-par,-par2)),as.integer(family-30),as.double(rep(0,length(u1))),PACKAGE='VineCopula')[[6]]
+ }else if(family %in% c(104,114,124,134,204,214,224,234)){ # Kan später ev. mal duch C-Code ersetzt werden
+ ## Hilfsterme für die Ableitung ###
+ ta<-function(t,par,par2,par3) {(par2*t)^par+(par3*(1-t))^par}
+ ######## Pickands A
+ A<-function(t,par,par2,par3) {(1-par3)*(1-t)+(1-par2)*t+ta(t,par,par2,par3)^(1/par)}
+
+ w<-function(u1,u2) {log(u2)/log(u1*u2)}
+ C<-function(u,v,par,par2,par3) {(u1*u2)^A(w(u1,u2),par,par2,par3)}
+
+ if(family==104)
+ {
+ par3=1
+ res=C(u1,u2,par,par2,par3)
+ }
+ else if(family==114)
+ {
+ par3=1
+ res=u1+u2-1+C(1-u1,1-u2,par,par2,par3)
+ }
+ else if(family==124)
+ {
+ par3=1
+ res=u2-C(1-u1,u2,-par,par2,par3)
+ }
+ else if(family==134)
+ {
+ par3=1
+ res=u1-C(u1,1-u2,-par,par2,par3)
+ }
+ else if(family==204)
+ {
+ par3=par2
+ par2=1
+ res=C(u1,u2,par,par2,par3)
+ }
+ else if(family==214)
+ {
+ par3=par2
+ par2=1
+ res=u1+u2-1+C(1-u1,1-u2,par,par2,par3)
+ }
+ else if(family==224)
+ {
+ par3=par2
+ par2=1
+ res=u2-C(1-u1,u2,-par,par2,par3)
+ }
+ else if(family==234)
+ {
+ par3=par2
+ par2=1
+ res=u1-C(u1,1-u2,-par,par2,par3)
+ }
}
return(res)
Modified: pkg/R/BiCopEst.r
===================================================================
--- pkg/R/BiCopEst.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopEst.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -22,7 +22,7 @@
if(length(u1)<2) stop("Number of observations has to be at least 2.")
if(any(u1>1) || any(u1<0)) stop("Data has be in the interval [0,1].")
if(any(u2>1) || any(u2<0)) stop("Data has be in the interval [0,1].")
- 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,51,61,71))) stop("Copula family not implemented.")
+ 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,51,61,71,104,114,124,134,204,214,224,234))) stop("Copula family not implemented.")
if(max.df<=2) stop("The upper bound for the degrees of freedom parameter has to be larger than 2.")
if(!is.list(max.BB)) stop("'max.BB' has to be a list.")
@@ -40,7 +40,7 @@
if(is.logical(se)==FALSE) stop("'se' has to be a logical variable (TRUE or FALSE).")
- if(method=="itau" && family %in% c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40))
+ if(method=="itau" && family %in% c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,104,114,124,134,204,214,224,234))
{
message("For two parameter copulas the estimation method 'itau' cannot be used. The method is automatically set to 'mle'.")
method="mle"
@@ -167,7 +167,7 @@
theta1=0
delta=0
- if(!(family%in%c(2,6,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40)))
+ if(!(family%in%c(2,6,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,104,114,124,134,204,214,224,234)))
{
theta1=theta
}
@@ -289,14 +289,77 @@
theta1=max(-1.5,-max((max.BB$BB8[1]+1.001)/2,1.001))
}
}
+ else if(family==104 || family==114)
+ {
+ if(tau<0)
+ {
+ print("The Tawn or survival Tawn copula cannot be used for negatively dependent data.")
+ delta=1
+ theta1=1.001
+ }
+ else
+ {
+ delta=0.5 # psi1
+ theta1=2
+ }
+ }
+ else if(family==124 || family==134)
+ {
+ if(tau>0)
+ {
+ print("The rotated Tawn copula cannot be used for positively dependent data.")
+ delta=1
+ theta1=-1.001
+ }
+ else
+ {
+ delta=0.5 # psi1
+ theta1=-2
+ }
+ }
+ else if(family==204 || family==214)
+ {
+ if(tau<0)
+ {
+ print("The Tawn2 or survival Tawn2 copula cannot be used for negatively dependent data.")
+ delta=1
+ theta1=1.001
+ }
+ else
+ {
+ delta=0.5 # psi1
+ theta1=2
+ }
+ }
+ else if(family==224 || family==234)
+ {
+ if(tau>0)
+ {
+ print("The rotated Tawn2 copula cannot be used for positively dependent data.")
+ delta=1
+ theta1=-1.001
+ }
+ else
+ {
+ delta=0.5 # psi1
+ theta1=-2
+ }
+ }
- if(family!=0)
+ if(family!=0 && family<100)
{
out=MLE_intern(cbind(u1,u2),c(theta1, delta),family=family,se,max.df,max.BB,weights)
theta=out$par
if(se==TRUE)
se1=out$se
}
+ else if(family!=0 && family>100) # New
+ {
+ out=MLE_intern_Tawn(cbind(u1,u2),c(theta1, delta),family=family,se)
+ theta=out$par
+ if(se==TRUE)
+ se1=out$se
+ }
}
@@ -421,11 +484,11 @@
if(is.null(weights))
{
- ll = .C("LL_mod",as.integer(family),as.integer(n),as.double(data[,2]),as.double(data[,1]),as.double(param[1]),as.double(param[2]),as.double(0),PACKAGE='VineCopula')[[7]]
+ ll = .C("LL_mod2",as.integer(family),as.integer(n),as.double(data[,1]),as.double(data[,2]),as.double(param[1]),as.double(param[2]),as.double(0),PACKAGE='VineCopula')[[7]]
}
else
{
- ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,2]), as.double(data[,1]), as.double(param[1]),as.double(param[2]), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
+ ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,1]), as.double(data[,2]), as.double(param[1]),as.double(param[2]), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
}
if(is.infinite(ll) || is.na(ll) || ll< -10^300) ll = -10^300
@@ -446,7 +509,7 @@
low = c(1.001,0.001)
up = max.BB$BB8
}else if(family == 27 | family==37){
- up = c(-1.001,-0.001)
+ up = c(-0.001,-1.001)
low = -max.BB$BB1
}else if(family == 28 | family==38){
up = c(-1.001,-1.001)
@@ -480,11 +543,11 @@
{
if(is.null(weights))
{
- ll = .C("LL_mod",as.integer(family),as.integer(n),as.double(data[,2]),as.double(data[,1]),as.double(param[1]),as.double(param[2]),as.double(0),PACKAGE='VineCopula')[[7]]
+ ll = .C("LL_mod2",as.integer(family),as.integer(n),as.double(data[,1]),as.double(data[,2]),as.double(param[1]),as.double(param[2]),as.double(0),PACKAGE='VineCopula')[[7]]
}
else
{
- ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,2]), as.double(data[,1]), as.double(param[1]),as.double(param[2]), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
+ ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,1]), as.double(data[,2]), as.double(param[1]),as.double(param[2]), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
}
if(is.infinite(ll) || is.na(ll) || ll< -10^10) ll = -10^10
@@ -523,11 +586,11 @@
{
if(is.null(weights))
{
- ll = .C("LL_mod",as.integer(family),as.integer(n),as.double(data[,2]),as.double(data[,1]),as.double(start.parm[1]),as.double(param[1]),as.double(0),PACKAGE='VineCopula')[[7]]
+ ll = .C("LL_mod2",as.integer(family),as.integer(n),as.double(data[,1]),as.double(data[,2]),as.double(start.parm[1]),as.double(param[1]),as.double(0),PACKAGE='VineCopula')[[7]]
}
else
{
- ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,2]), as.double(data[,1]), as.double(start.parm[1]),as.double(param[1]), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
+ ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,1]), as.double(data[,2]), as.double(start.parm[1]),as.double(param[1]), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
}
if(is.infinite(ll) || is.na(ll) || ll< -10^300) ll = -10^300
@@ -563,11 +626,11 @@
{
if(is.null(weights))
{
- ll = .C("LL_mod",as.integer(family),as.integer(n),as.double(data[,2]), as.double(data[,1]), as.double(param),as.double(0), as.double(0),PACKAGE='VineCopula')[[7]]
+ ll = .C("LL_mod2",as.integer(family),as.integer(n),as.double(data[,1]), as.double(data[,2]), as.double(param),as.double(0), as.double(0),PACKAGE='VineCopula')[[7]]
}
else
{
- ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,2]), as.double(data[,1]), as.double(param[1]),as.double(0), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
+ ll = .C("LL_mod_seperate",as.integer(family),as.integer(n),as.double(data[,1]), as.double(data[,2]), as.double(param[1]),as.double(0), as.double(rep(0,n)),PACKAGE='VineCopula')[[7]]%*%weights
}
if(is.infinite(ll) || is.na(ll) || ll< -10^300) ll = -10^300
@@ -707,6 +770,69 @@
}
+# New for Tawn
+
+MLE_intern_Tawn <-
+function(data,start.parm,family,se=FALSE)
+{
+
+ n = dim(data)[1]
+ tau <- fasttau(data[,1],data[,2])
+
+ if(family==104 || family==114 || family==204 || family==214)
+ {
+ parlower<-c(1.001,max(tau,0.0001))
+ parupper<-c(20,min(tau+0.1,0.99))
+ }
+ else if(family==124 || family==134 || family==224 || family==234)
+ {
+ parlower<-c(-20,max(-tau,0.0001))
+ parupper<-c(-1.001,min(-tau+0.1,0.99))
+ }
+
+ # Hier fehlt noch die log-likelihood Funktion
+ loglikfunc = function(param)
+ {
+ ll = .C("LL_mod2",as.integer(family),as.integer(n),as.double(data[,1]), as.double(data[,2]), as.double(param[1]),as.double(param[2]), as.double(0),PACKAGE='VineCopula')[[7]]
+ if(is.infinite(ll) || is.na(ll) || ll< -10^300) ll = -10^300
+ #print(param)
+ #print(ll)
+ return(ll)
+ }
+
+ out=list()
+ #print(start.parm)
+ if(se == TRUE)
+ {
+ optimout=optim(par=start.parm,fn=loglikfunc,method=c("L-BFGS-B"),lower=parlower,upper=parupper,control=list(fnscale=-1,maxit=500), hessian=TRUE)
+ if(det(optimout$hessian)==0){
+ var = diag(1,dim(optimout$hessian)[1])
+ }else{
+ var = (-solve(optimout$hessian))
+ }
+
+ out$se = sqrt(diag(var))
+ }
+ else
+ {
+ optimout=optim(par=start.parm,fn=loglikfunc,method=c("L-BFGS-B"),lower=parlower,upper=parupper,control=list(fnscale=-1,maxit=500))
+ }
+
+ out$par = optimout$par
+ out$value=optimout$value
+ return(out)
+}
+
+
+
+
+
+
+
+
+
+
+
fasttau<- function(x, y,weights=NA)
{
if(any(is.na(weights))){
Modified: pkg/R/BiCopGofTest.r
===================================================================
--- pkg/R/BiCopGofTest.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopGofTest.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -1,5 +1,8 @@
-BiCopGofTest<-function(u1, u2, family, par=0, par2=0, method="white", max.df=30, B=100, level=0.05)
+BiCopGofTest<-function(u1, u2, family, par=0, par2=0, method="white", max.df=30, B=100)
{
+ if(method=="White") method="white"
+ if(method=="Kendall") method="kendall"
+
if(is.null(u1)==TRUE || is.null(u2)==TRUE) stop("u1 and/or u2 are not set or have length zero.")
if(length(u1)!=length(u2)) stop("Lengths of 'u1' and 'u2' do not match.")
if(any(u1>1) || any(u1<0)) stop("Data has be in the interval [0,1].")
@@ -38,7 +41,7 @@
}
if(family==2 && method=="kendall") stop("The goodness-of-fit test based on Kendall's process is not implemented for the t-copula.")
if(family%in%c(7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40) && method=="white") stop("The goodness-of-fit test based on White's information matrix equality is not implemented for the BB copulas.")
- if((level < 0 || level > 1) && method=="kendall") stop("Significance level has to be between 0 and 1.")
+ #if((level < 0 || level > 1) && method=="kendall") stop("Significance level has to be between 0 and 1.")
T=length(u1)
@@ -123,6 +126,75 @@
}
out=list(p.value=pvalue,statistic=test)
}
+ else if(method=="IR")
+ {
+ # Information ratio GOF
+ # Step 1: maximum likelihood estimation
+ if(family==2)
+ {
+ if(par==0)
+ {
+ pars=BiCopEst(u1,u2,family=family,method="mle",max.df=max.df)
+ theta=pars$par
+ nu=pars$par2
+ print(theta)
+ print(nu)
+ }
+ else
+ {
+ theta=par
+ nu=par2
+ }
+ }
+ else
+ {
+ nu=0
+ theta=BiCopEst(u1,u2,family=family,method="mle")$par
+ }
+
+ # Step 2: Calculation of Hesse and gradient
+ if(family==2)
+ {
+ grad=c(0,0)
+ rho_teil=f_rho(u1,u2,theta,nu)
+ nu_teil=f_nu(u1,u2,theta,nu)
+ rho_nu_teil=f_rho_nu(u1,u2,theta,nu)
+ H = matrix(c(rho_teil,rho_nu_teil,rho_nu_teil,nu_teil),2,2) # Hesse matrix
+ grad[1]=BiCopDeriv(u1,u2,family=family,par=theta,par2=nu,deriv="par", log=TRUE)
+ grad[2]=BiCopDeriv(u1,u2,family=family,par=theta,par2=nu,deriv="par2", log=TRUE)
+ C=grad%*%t(grad)
+ }
+ else
+ {
+ d=rep(0,T)
+ for(t in 1:T)
+ {
+ b=BiCopPDF(u1[t],u2[t],family,theta,nu)
+ d[t]=BiCopDeriv2(u1[t],u2[t],family=family,par=theta,par2=nu, deriv="par")/b
+ }
+ H=mean(d)
+ C=BiCopDeriv(u1,u2,family=family,par=theta,par2=nu,deriv="par", log=TRUE)
+ }
+ Phi=-solve(H)%*%C
+ IR=trace(Phi)/dim(H)[1] #Zwischenergebnis
+
+ #Bootstrap procedure
+ if(B==0)
+ {
+ out=list(IR=IR, p.value=NULL)
+ }
+ else
+ {
+ IR_boot=boot.IR(family,theta,nu,B,length(u1))
+ sigma2=var(IR_boot)
+ IR_new=((IR-1)/sqrt(sigma2))^2
+ IR_boot=((IR_boot-1)/sqrt(sigma2))^2
+ p.value = mean(IR_boot>=IR_new)
+
+ out=list(IR=IR,p.value=p.value)
+ }
+
+ }
else if(method=="kendall")
{
if(family %in% c(13,14,16,17,18,19,20)){
@@ -161,9 +233,9 @@
sn.obs<-ostat$Sn
tn.obs<-ostat$Tn
- k<-as.integer((1-level)*B)
- sn.critical <- sn.boot[k] # critical value of test at level 0.05
- tn.critical <- tn.boot[k] # critical value of test at level 0.05
+ #k<-as.integer((1-level)*B)
+ #sn.critical <- sn.boot[k] # critical value of test at level 0.05
+ #tn.critical <- tn.boot[k] # critical value of test at level 0.05
pv.sn<-sapply(sn.obs,function(x) (1/B)*length(which(sn.boot[1:B]>=x))) # P-value of Sn
@@ -251,7 +323,7 @@
sam.par<-suppressWarnings({BiCopEst(sam[,1], sam[,2],family=fam)}) # parameter estimation of sample data
sim<-BiCopSim(10000,fam,sam.par$par,sam.par$par2) # generate data for the simulation of theo. K(t)
- #par2 muss auf einen Integer gesetzt werden für mvtnorm
+ #par2 muss auf einen Integer gesetzt werden f?r mvtnorm
param$par2=round(param$par2)
cormat = matrix(c(1,param$par,param$par,1),2,2)
@@ -436,3 +508,50 @@
return(out)
}
+
+############################
+
+# bootstrap for IR
+
+boot.IR<-function(family,theta,nu,B,n)
+{
+ #theta und nu sind die geschaetzten Parameter
+ IR=rep(0,B)
+ for(i in 1:B)
+ {
+ sam=BiCopSim(n,family,theta,nu)
+ sam.par<-BiCopEst(sam[,1], sam[,2],family=family) # parameter estimation of sample data
+ if(family==2)
+ {
+ theta2=sam.par[1]
+ nu2=sam.par[2]
+ grad=c(0,0)
+ rho_teil=f_rho(sam[,1],sam[,2],theta2,nu2)
+ nu_teil=f_nu(sam[,1],sam[,2],theta2,nu2)
+ rho_nu_teil=f_rho_nu(sam[,1],sam[,2],theta2,nu2)
+ H = matrix(c(rho_teil,rho_nu_teil,rho_nu_teil,nu_teil),2,2) # Hesse matrix
+ grad[1]=BiCopDeriv(sam[,1],sam[,2],family=family,par=theta2,par2=nu2,deriv="par", log=TRUE)
+ grad[2]=BiCopDeriv(sam[,1],sam[,2],family=family,par=theta2,par2=nu2,deriv="par2", log=TRUE)
+ C=grad%*%t(grad)
+ }
+ else
+ {
+ theta2=sam.par
+ nu2=0
+ d=rep(0,T)
+ for(t in 1:T)
+ {
+ b=BiCopPDF(sam[t,1],sam[t,2],family,theta,nu)
+ d[t]=BiCopDeriv2(sam[t,1],sam[t,2],family=family,par=theta,par2=nu, deriv="par")/b
+ }
+ H=mean(d)
+ C=BiCopDeriv(sam[,1],sam[,2],family=family,par=theta2,par2=nu2,deriv="par", log=TRUE)
+ }
+ Phi=-solve(H)%*%C
+ IR[i]=trace(Phi)/dim(H)[1]
+ }
+
+ return(IR)
+}
+
+
Modified: pkg/R/BiCopHfunc.r
===================================================================
--- pkg/R/BiCopHfunc.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopHfunc.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -17,10 +17,10 @@
if(is.null(u1)==TRUE || is.null(u2)==TRUE) stop("u1 and/or u2 are not set or have length zero.")
if(length(u1) != length(u2)) stop("Lengths of 'u1' and 'u2' do not match.")
if(any(u1>1) || any(u1<0)) stop("Data has be in the interval [0,1].")
- if(any(u2>1) || any(u2<0)) stop("Data has be in the interval [0,1].")
- 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,43,44, 41,51,61,71))) stop("Copula family not implemented.")
- if(c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40) %in% family && par2==0) stop("For t-, BB1, BB6, BB7 and BB8 copulas, 'par2' must be set.")
- if(c(1,3,4,5,6,13,14,16,23,24,26,33,34,36,41,51,61,71) %in% family && length(par)<1) stop("'par' not set.")
+ if(any(u2>1) || any(u2<0)) stop("Data has be in the interval [0,1].")
+ 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(2,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 t-, 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.")
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==2 && par2<=2) stop("The degrees of freedom parameter of the t-copula has to be larger than 2.")
@@ -49,6 +49,10 @@
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==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].")
n=length(u1)
Modified: pkg/R/BiCopMetaContour.r
===================================================================
--- pkg/R/BiCopMetaContour.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopMetaContour.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -213,6 +213,7 @@
}
+
#################################################################################
# cop.pdf #
# Input: #
@@ -462,6 +463,62 @@
d2=1-u2
return( pmax(a * d2 * (((12 - 9 * d1) * d1 - 3) * d2 + d1 * (6 * d1 - 8) + 2) + b * (d2 * ((d1 * (9 * d1 - 12) + 3) * d2 + (12 - 6 * d1) * d1 - 4) - 2 * d1 + 1) + 1,0) )
}
+ else if(copula==104) # Tawn copula (psi2 fix)
+ {
+ par=param[1]
+ par2=param[2]
+ fam=104
+ return(BiCopPDF(u1,u2,fam,par,par2))
+ }
+ else if(copula==114) # Tawn copula (psi2 fix)
+ {
+ par=param[1]
+ par2=param[2]
+ fam=104
+ return(BiCopPDF(1-u1,1-u2,fam,par,par2))
+ }
+ else if(copula==124) # Tawn copula (psi2 fix)
+ {
+ par=-param[1]
+ par2=param[2]
+ fam=104
+ return(BiCopPDF(1-u1,u2,fam,par,par2))
+ }
+ else if(copula==134) # Tawn copula (psi2 fix)
+ {
+ par=-param[1]
+ par2=param[2]
+ fam=104
+ return(BiCopPDF(u1,1-u2,fam,par,par2))
+ }
+ else if(copula==204) # Tawn copula (psi1 fix)
+ {
+ par=param[1]
+ par2=param[2]
+ fam=204
+ return(BiCopPDF(u1,u2,fam,par,par2))
+ }
+ else if(copula==214) # Tawn copula (psi1 fix)
+ {
+ par=param[1]
+ par2=param[2]
+ fam=204
+ return(BiCopPDF(1-u1,1-u2,fam,par,par2))
+ }
+ else if(copula==224) # Tawn copula (psi1 fix)
+ {
+ par=-param[1]
+ par2=param[2]
+ fam=204
+ return(BiCopPDF(1-u1,u2,fam,par,par2))
+ }
+ else if(copula==234) # Tawn copula (psi1 fix)
+ {
+ par=-param[1]
+ par2=param[2]
+ fam=204
+ return(BiCopPDF(u1,1-u2,fam,par,par2))
+ }
}
@@ -514,8 +571,8 @@
if(is.null(u1)==FALSE && (any(u1>1) || any(u1<0))) stop("Data has be in the interval [0,1].")
if(is.null(u2)==FALSE && (any(u2>1) || any(u2<0))) stop("Data has be in the interval [0,1].")
#if(length(u1)!=length(u2)) stop("Lengths of 'u1' and 'u2' do not match.")
- 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, "emp"))) stop("Copula family not implemented.")
- if(c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,42,52,62,72) %in% family && par2==0) stop("For t-, BB1 and BB7 copulas, 'par2' must be set.")
+ 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, "emp"))) stop("Copula family not implemented.")
+ if(c(2,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 t-, 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.")
# size sollte nicht zu gross sein
@@ -562,6 +619,10 @@
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==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(PLOT!=TRUE && PLOT!=FALSE) stop("The parameter 'PLOT' has to be set to 'TRUE' or 'FALSE'.")
@@ -616,7 +677,7 @@
if(family!="emp")
{
- if(family %in% c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,42,52,62,72))
+ if(family %in% c(2,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))
z <- matrix(data=meta.dens(x1=rep(x=x, each=size), x2=rep(x=y, times=size), param=c(par,par2), copula=family,
margins=margins, margins.par=margins.par), nrow=size, byrow=TRUE)
else
Modified: pkg/R/BiCopName.r
===================================================================
--- pkg/R/BiCopName.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopName.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -43,6 +43,14 @@
else if(family==51) fam="1-par AS180"
else if(family==61) fam="1-par AS90"
else if(family==71) fam="1-par AS270"
+ else if(family==104) fam="Tawn"
+ else if(family==114) fam="Tawn180"
+ else if(family==124) fam="Tawn90"
+ else if(family==134) fam="Tawn270"
+ else if(family==204) fam="Tawn2"
+ else if(family==214) fam="Tawn2_180"
+ else if(family==224) fam="Tawn2_90"
+ else if(family==234) fam="Tawn2_270"
else stop("Family not implemented.")
}
else # langer Name
@@ -84,6 +92,14 @@
else if(family==51) fam="Rotated 1-parametric asymmetric 180 degree"
else if(family==61) fam="Rotated 1-parametric asymmetric 90 degree"
else if(family==71) fam="Rotated 1-parametric asymmetric 270 degree"
+ else if(family==104) fam="Tawn"
+ else if(family==114) fam="Rotated Tawn 180 degrees"
+ else if(family==124) fam="Rotated Tawn 90 degrees"
+ else if(family==134) fam="Rotated Tawn 270 degrees"
+ else if(family==204) fam="Tawn2"
+ else if(family==214) fam="Rotated Tawn2 180 degrees"
+ else if(family==224) fam="Rotated Tawn2 90 degrees"
+ else if(family==234) fam="Rotated Tawn2 270 degrees"
else stop("Family not implemented.")
}
}
@@ -125,6 +141,14 @@
else if(family=="1-par AS180" || family=="Rotated 1-parametric asymmetric 180 degree") fam=51
else if(family=="1-par AS90" || family=="Rotated 1-parametric asymmetric 90 degree") fam=61
else if(family=="1-par AS270" || family=="Rotated 1-parametric asymmetric 270 degree") fam=71
+ else if(family=="Tawn") fam=104
+ else if(family=="Tawn180" || family=="Rotated Tawn 180 degrees") fam=114
+ else if(family=="Tawn90" || family=="Rotated Tawn 90 degrees") fam=124
+ else if(family=="Tawn270" || family=="Rotated Tawn 270 degrees") fam=134
+ else if(family=="Tawn2") fam=204
+ else if(family=="Tawn2_180" || family=="Rotated Tawn2 180 degrees") fam=214
+ else if(family=="Tawn2_90" || family=="Rotated Tawn2 90 degrees") fam=224
+ else if(family=="Tawn2_270" || family=="Rotated Tawn2 270 degrees") fam=234
else stop("Family not implemented.")
}
Modified: pkg/R/BiCopPDF.r
===================================================================
--- pkg/R/BiCopPDF.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopPDF.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -4,8 +4,8 @@
if(length(u1)!=length(u2)) stop("Lengths of 'u1' and 'u2' do not match.")
if(any(u1>1) || any(u1<0)) stop("Data has be in the interval [0,1].")
if(any(u2>1) || any(u2<0)) stop("Data has be in the interval [0,1].")
- 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,51,61,71))) stop("Copula family not implemented.")
- if(family %in% c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40) && par2==0) stop("For t-, BB1, BB6, BB7 and BB8 copulas, 'par2' must be set.")
+ 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,51,61,71,104,114,124,134,204,214,224,234))) stop("Copula family not implemented.")
+ if(family %in% c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40,104,114,124,134,204,214,224,234) && par2==0) stop("For t-, BB1, BB6, BB7, BB8 and Tawn copulas, 'par2' must be set.")
if(family %in% c(1,3,4,5,6,13,14,16,23,24,26,33,34,36,41,51,61,71) && length(par)<1) stop("'par' not 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).")
@@ -35,6 +35,10 @@
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==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].")
coplik = .C("LL_mod_seperate",as.integer(family),as.integer(length(u1)),as.double(u1),as.double(u2),as.double(par),as.double(par2),as.double(rep(0,length(u1))),PACKAGE='VineCopula')[[7]]
Added: pkg/R/BiCopPar2Beta.r
===================================================================
--- pkg/R/BiCopPar2Beta.r (rev 0)
+++ pkg/R/BiCopPar2Beta.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -0,0 +1,6 @@
+BiCopPar2Beta <- function(family,par,par2=0)
+{
+ blomBeta=4*BiCopCDF(0.5,0.5,family,par,par2)-1
+
+ return(blomBeta)
+}
\ No newline at end of file
Modified: pkg/R/BiCopPar2TailDep.r
===================================================================
--- pkg/R/BiCopPar2TailDep.r 2013-09-18 13:08:33 UTC (rev 24)
+++ pkg/R/BiCopPar2TailDep.r 2013-10-09 13:49:16 UTC (rev 25)
@@ -1,8 +1,8 @@
BiCopPar2TailDep<-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))) stop("Copula family not implemented.")
- if(c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40) %in% family && par2==0) stop("For t-, BB1, BB6, BB7 and BB8 copulas, 'par2' must be set.")
- if(c(1,3,4,5,6,13,14,16,23,24,26,33,34,36) %in% family && length(par)<1) stop("'par' not set.")
+ 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(2,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 t-, 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.")
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).")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 25
Mehr Informationen über die Mailingliste Vinecopula-commits