[Vinecopula-commits] r6 - / pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Do Apr 4 19:33:11 CEST 2013
Author: ben_graeler
Date: 2013-04-04 19:33:10 +0200 (Thu, 04 Apr 2013)
New Revision: 6
Modified:
/
pkg/R/RVineAIC.r
pkg/R/RVineGrad.r
pkg/R/RVineHessian.r
pkg/R/RVineLogLik.r
pkg/R/RVineMLE.R
pkg/R/RVinePar2Tau.r
pkg/R/RVineSeqEst.R
pkg/R/RVineSim.R
pkg/R/RVineTreePlot.r
Log:
- replaced is(RVM) != "RVineMatrix" with !is(RVM, "RVineMatrix") in several if-statetments to avoid warnings in case "RVM" possesses several S3-classes.
Property changes on:
___________________________________________________________________
Added: svn:ignore
+ .Rproj.user
.Rhistory
.RData
Modified: pkg/R/RVineAIC.r
===================================================================
--- pkg/R/RVineAIC.r 2013-04-04 09:24:02 UTC (rev 5)
+++ pkg/R/RVineAIC.r 2013-04-04 17:33:10 UTC (rev 6)
@@ -11,7 +11,7 @@
n<-d
N<-T
if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
- if(is(RVM) != "RVineMatrix") stop("'RVM' has to be an RVineMatrix object.")
+ if(!is(RVM,"RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
npar = sum(RVM$family >= 1, na.rm=TRUE) + sum(RVM$family %in% c(2,7:10,17:20,27:30,37:40),na.rm=TRUE)
npar_pair = (RVM$family>=1)+(RVM$family%in%c(2,7:10,17:20,27:30,37:40))
Property changes on: pkg/R/RVineAIC.r
___________________________________________________________________
Added: svn:eol-style
+ LF
Modified: pkg/R/RVineGrad.r
===================================================================
--- pkg/R/RVineGrad.r 2013-04-04 09:24:02 UTC (rev 5)
+++ pkg/R/RVineGrad.r 2013-04-04 17:33:10 UTC (rev 6)
@@ -1,150 +1,150 @@
-#################################################################
-# #
-# RVineGrad #
-# #
-# Function to calculate the derivative of one #
-# pair-copula in an R-vine #
-# #
-# Input: #
-# data data set #
-# RVM R-Vine matrix object #
-# calcupdate array of Update-Matrices (output of RVineMatrixUpdate) #
-# par, par2 Copula parameter stored in an RVM-Matrix #
-# start.V log-liklihoods (output of RVineLogLik) #
-# #
-# Output: #
-# gradient gradient of the R-vine #
-#################################################################
-
-RVineGrad <-function(data,RVM,par=RVM$par,par2=RVM$par2,start.V=NA, posParams=(RVM$family > 0))
-{
-
-if(any(!(RVM$family %in% c(0,1:6,13,14,16,23,24,26,33,34,36)))) stop("Copula family not implemented.")
-
-if(is.vector(data)){
- data = t(as.matrix(data))
- }else{
- data=as.matrix(data)
- }
-
- if(any(data>1) || any(data<0)) stop("Data has be in the interval [0,1].")
- d=dim(data)[2]
- T=dim(data)[1]
- n<-d
- N<-T
- if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
- if(is(RVM) != "RVineMatrix") stop("'RVM' has to be an RVineMatrix object.")
-
-
-#if(any(is.na(calcupdate)))
-#{
-# n=dim(RVM)
-# calcupdate=array(0,dim=c(n,n,n,n))
-# for(i in (n-1):1){
-# for(k in n:(i+1)){
-# calcupdate[, ,k,i ]=RVineMatrixUpdate(RVM,k,i)
-# }
-# }
-#}
-
-o = diag(RVM$Matrix)
- if(any(o != length(o):1))
- {
- oldRVM = RVM
- RVM = normalizeRVineMatrix(RVM)
- #RVM = getFromNamespace("normalizeRVineMatrix","VineCopula")(RVM)
- data = data[,o[length(o):1]]
- }
-
- if(any(is.na(start.V)))
- {
- loglik=RVineLogLik(data,RVM,par=par, par2=par2,separate=TRUE)
- V=loglik$V
-
- }
- else
- {
- V = start.V
- V$value[V$value %in% c(NA,NaN,-Inf)]<- -1e10
- if(any(is.na(V$value))) message("NA in LogL call")
- }
-
-
- ll=as.vector(V$value)
- vv=as.vector(V$direct)
- vv2=as.vector(V$indirect)
- #calcup=as.vector(calcupdate)
-
- w1=as.vector(RVM$family)
- w1[is.na(w1)]<-0
- th=as.vector(par)
- th[is.na(th)]<-0
- th2=as.vector(par2)
- th2[is.na(th2)]<-0
- condirect=as.vector(as.numeric(RVM$CondDistr$direct))
- conindirect=as.vector(as.numeric(RVM$CondDistr$indirect))
- maxmat=as.vector(RVM$MaxMat)
- matri=as.vector(RVM$Matrix)
- matri[is.na(matri)]<-0
- maxmat[is.na(maxmat)]<-0
- condirect[is.na(condirect)]<-0
- conindirect[is.na(conindirect)]<-0
- # tilde_vdirect_array=array(0,dim=c(n,n,N,n,n))
- # tilde_vindirect_array=array(0,dim=c(n,n,N,n,n))
- # tilde_value_array=array(0,dim=c(n,n,N,n,n))
-
-
- out=rep(0,sum(posParams[lower.tri(posParams, diag=FALSE)])+sum(w1==2))
-
-
-out <- .C("VineLogLikRvineGradient",
- as.integer(T),
- as.integer(d),
- as.integer(w1),
- as.integer(maxmat),
- as.integer(matri),
- as.integer(condirect),
- as.integer(conindirect),
- as.double(th),
- as.double(th2),
- as.double(data),
- as.double(out),
- as.double(ll),
- as.double(vv),
- as.double(vv2),
- #as.integer(calcup),
- as.integer(as.vector(posParams)),
- #as.double(as.vector(tilde_vdirect_array)),
- #as.double(as.vector(tilde_vindirect_array)),
- #as.double(as.vector(tilde_value_array)),
- PACKAGE = 'VineCopula')
-
-
- gradient2 <- out[[11]]
- gradient2[gradient2 %in% c(NA,NaN,-Inf)] <- -1e10
-
- dd=sum(RVM$family>0)
- tt=sum(w1==2)
- grad1=gradient2[1:dd]
- gradient=grad1[dd:1]
- if(tt>0)
- {
- grad2=gradient2[(dd+1):(dd+tt)]
- gradient=c(gradient,grad2[tt:1])
- }
-
-
- #tilde_vdirect=out[[16]]
- #tilde_vindirect=out[[17]]
- #tilde_value=out[[18]]
- #V$tilde_direct = array(tilde_vdirect,dim=c(n,n,N,n,n))
- #V$tilde_indirect = array(tilde_vindirect,dim=c(n,n,N,n,n))
- #V$tilde_value = array(tilde_value,dim=c(n,n,N,n,n))
-
-
-
-
-#out2=list(gradient=gradient,V=V)
-out2=list(gradient=gradient)
-return(out2)
-}
+#################################################################
+# #
+# RVineGrad #
+# #
+# Function to calculate the derivative of one #
+# pair-copula in an R-vine #
+# #
+# Input: #
+# data data set #
+# RVM R-Vine matrix object #
+# calcupdate array of Update-Matrices (output of RVineMatrixUpdate) #
+# par, par2 Copula parameter stored in an RVM-Matrix #
+# start.V log-liklihoods (output of RVineLogLik) #
+# #
+# Output: #
+# gradient gradient of the R-vine #
+#################################################################
+
+RVineGrad <-function(data,RVM,par=RVM$par,par2=RVM$par2,start.V=NA, posParams=(RVM$family > 0))
+{
+
+if(any(!(RVM$family %in% c(0,1:6,13,14,16,23,24,26,33,34,36)))) stop("Copula family not implemented.")
+
+if(is.vector(data)){
+ data = t(as.matrix(data))
+ }else{
+ data=as.matrix(data)
+ }
+
+ if(any(data>1) || any(data<0)) stop("Data has be in the interval [0,1].")
+ d=dim(data)[2]
+ T=dim(data)[1]
+ n<-d
+ N<-T
+ if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
+ if(!is(RVM,"RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
+
+
+#if(any(is.na(calcupdate)))
+#{
+# n=dim(RVM)
+# calcupdate=array(0,dim=c(n,n,n,n))
+# for(i in (n-1):1){
+# for(k in n:(i+1)){
+# calcupdate[, ,k,i ]=RVineMatrixUpdate(RVM,k,i)
+# }
+# }
+#}
+
+o = diag(RVM$Matrix)
+ if(any(o != length(o):1))
+ {
+ oldRVM = RVM
+ RVM = normalizeRVineMatrix(RVM)
+ #RVM = getFromNamespace("normalizeRVineMatrix","VineCopula")(RVM)
+ data = data[,o[length(o):1]]
+ }
+
+ if(any(is.na(start.V)))
+ {
+ loglik=RVineLogLik(data,RVM,par=par, par2=par2,separate=TRUE)
+ V=loglik$V
+
+ }
+ else
+ {
+ V = start.V
+ V$value[V$value %in% c(NA,NaN,-Inf)]<- -1e10
+ if(any(is.na(V$value))) message("NA in LogL call")
+ }
+
+
+ ll=as.vector(V$value)
+ vv=as.vector(V$direct)
+ vv2=as.vector(V$indirect)
+ #calcup=as.vector(calcupdate)
+
+ w1=as.vector(RVM$family)
+ w1[is.na(w1)]<-0
+ th=as.vector(par)
+ th[is.na(th)]<-0
+ th2=as.vector(par2)
+ th2[is.na(th2)]<-0
+ condirect=as.vector(as.numeric(RVM$CondDistr$direct))
+ conindirect=as.vector(as.numeric(RVM$CondDistr$indirect))
+ maxmat=as.vector(RVM$MaxMat)
+ matri=as.vector(RVM$Matrix)
+ matri[is.na(matri)]<-0
+ maxmat[is.na(maxmat)]<-0
+ condirect[is.na(condirect)]<-0
+ conindirect[is.na(conindirect)]<-0
+ # tilde_vdirect_array=array(0,dim=c(n,n,N,n,n))
+ # tilde_vindirect_array=array(0,dim=c(n,n,N,n,n))
+ # tilde_value_array=array(0,dim=c(n,n,N,n,n))
+
+
+ out=rep(0,sum(posParams[lower.tri(posParams, diag=FALSE)])+sum(w1==2))
+
+
+out <- .C("VineLogLikRvineGradient",
+ as.integer(T),
+ as.integer(d),
+ as.integer(w1),
+ as.integer(maxmat),
+ as.integer(matri),
+ as.integer(condirect),
+ as.integer(conindirect),
+ as.double(th),
+ as.double(th2),
+ as.double(data),
+ as.double(out),
+ as.double(ll),
+ as.double(vv),
+ as.double(vv2),
+ #as.integer(calcup),
+ as.integer(as.vector(posParams)),
+ #as.double(as.vector(tilde_vdirect_array)),
+ #as.double(as.vector(tilde_vindirect_array)),
+ #as.double(as.vector(tilde_value_array)),
+ PACKAGE = 'VineCopula')
+
+
+ gradient2 <- out[[11]]
+ gradient2[gradient2 %in% c(NA,NaN,-Inf)] <- -1e10
+
+ dd=sum(RVM$family>0)
+ tt=sum(w1==2)
+ grad1=gradient2[1:dd]
+ gradient=grad1[dd:1]
+ if(tt>0)
+ {
+ grad2=gradient2[(dd+1):(dd+tt)]
+ gradient=c(gradient,grad2[tt:1])
+ }
+
+
+ #tilde_vdirect=out[[16]]
+ #tilde_vindirect=out[[17]]
+ #tilde_value=out[[18]]
+ #V$tilde_direct = array(tilde_vdirect,dim=c(n,n,N,n,n))
+ #V$tilde_indirect = array(tilde_vindirect,dim=c(n,n,N,n,n))
+ #V$tilde_value = array(tilde_value,dim=c(n,n,N,n,n))
+
+
+
+
+#out2=list(gradient=gradient,V=V)
+out2=list(gradient=gradient)
+return(out2)
+}
Property changes on: pkg/R/RVineGrad.r
___________________________________________________________________
Added: svn:eol-style
+ LF
Modified: pkg/R/RVineHessian.r
===================================================================
--- pkg/R/RVineHessian.r 2013-04-04 09:24:02 UTC (rev 5)
+++ pkg/R/RVineHessian.r 2013-04-04 17:33:10 UTC (rev 6)
@@ -1,91 +1,91 @@
-RVineHessian <-function(data,RVM)
-{
-
- if(any(!(RVM$family %in% c(0,1:6,13,14,16,23,24,26,33,34,36)))) stop("Copula family not implemented.")
-
- if(is.vector(data)){
- data = t(as.matrix(data))
- }else{
- data=as.matrix(data)
- }
- if(any(data>1) || any(data<0)) stop("Data has be in the interval [0,1].")
- if(is.null(dim(data)))
- {
- d=length(data)
- T=1
- }
- else
- {
- d=dim(data)[2]
- T=dim(data)[1]
- }
- n<-d
- N<-T
- if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
- if(is(RVM) != "RVineMatrix") stop("'RVM' has to be an RVineMatrix object.")
-
-
- o = diag(RVM$Matrix)
- if(any(o != length(o):1))
- {
- oldRVM = RVM
- #RVM = normalizeRVineMatrix(RVM)
- RVM = getFromNamespace("normalizeRVineMatrix","VineCopula")(RVM)
- data = data[,o[length(o):1]]
- }
-
- dd=d*(d-1)/2
- tt=sum(RVM$family==2)
- hessian=matrix(0,dd+tt,dd+tt)
- subhess=matrix(0,dd+tt,dd+tt)
- der=matrix(0,dd+tt,dd+tt)
- subder=matrix(0,dd+tt,dd+tt)
-
- out <- .C("hesse",
- as.integer(T),
- as.integer(d),
- as.integer(as.vector(RVM$family)),
- as.integer(as.vector(RVM$MaxMat)),
- as.integer(as.vector(RVM$Matrix)),
- as.integer(as.vector(RVM$CondDistr$direct)),
- as.integer(as.vector(RVM$CondDistr$indirect)),
- as.double(as.vector(RVM$par)),
- as.double(as.vector(RVM$par2)),
- as.double(as.vector(data)),
- as.double(as.vector(hessian)),
- as.double(as.vector(subhess)),
- as.double(as.vector(der)),
- as.double(as.vector(subder)),
- PACKAGE = 'VineCopula')
-
- hessian=matrix(out[[11]],dd+tt,dd+tt)
- subhess=matrix(out[[12]],dd+tt,dd+tt)
- der=matrix(out[[13]],dd+tt,dd+tt)
- subder=matrix(out[[14]],dd+tt,dd+tt)
-
- # der[1:dd,1:dd]=der[dd:1,dd:1]
- # if(tt>0)
- # {
- # der[(dd+1):(dd+tt),1:dd]=der[(dd+tt):(dd+1),dd:1]
- # der[1:dd,(dd+1):(dd+tt)]=der[dd:1,(dd+tt):(dd+1)]
- # der[(dd+1):(dd+tt),(dd+1):(dd+tt)]=der[(dd+tt):(dd+1),(dd+tt):(dd+1)]
- # }
- # hessian[1:dd,1:dd]=hessian[dd:1,dd:1]
- # if(tt>0)
- # {
- # hessian[(dd+1):(dd+tt),1:dd]=hessian[(dd+tt):(dd+1),dd:1]
- # hessian[1:dd,(dd+1):(dd+tt)]=hessian[dd:1,(dd+tt):(dd+1)]
- # hessian[(dd+1):(dd+tt),(dd+1):(dd+tt)]=hessian[(dd+tt):(dd+1),(dd+tt):(dd+1)]
- # }
-
- test=apply(hessian,2,function(x) max(abs(x)))
- hessian=hessian[test>0,test>0]
- subhess=subhess[test>0,test>0]
- der=der[test>0,test>0]
- subder=subder[test>0,test>0]
-
-
- out=list(hessian=hessian, der=der)
-
-return(out)
-}
+RVineHessian <-function(data,RVM)
+{
+
+ if(any(!(RVM$family %in% c(0,1:6,13,14,16,23,24,26,33,34,36)))) stop("Copula family not implemented.")
+
+ if(is.vector(data)){
+ data = t(as.matrix(data))
+ }else{
+ data=as.matrix(data)
+ }
+ if(any(data>1) || any(data<0)) stop("Data has be in the interval [0,1].")
+ if(is.null(dim(data)))
+ {
+ d=length(data)
+ T=1
+ }
+ else
+ {
+ d=dim(data)[2]
+ T=dim(data)[1]
+ }
+ n<-d
+ N<-T
+ if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
+ if(!is(RVM,"RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
+
+
+ o = diag(RVM$Matrix)
+ if(any(o != length(o):1))
+ {
+ oldRVM = RVM
+ #RVM = normalizeRVineMatrix(RVM)
+ RVM = getFromNamespace("normalizeRVineMatrix","VineCopula")(RVM)
+ data = data[,o[length(o):1]]
+ }
+
+ dd=d*(d-1)/2
+ tt=sum(RVM$family==2)
+ hessian=matrix(0,dd+tt,dd+tt)
+ subhess=matrix(0,dd+tt,dd+tt)
+ der=matrix(0,dd+tt,dd+tt)
+ subder=matrix(0,dd+tt,dd+tt)
+
+ out <- .C("hesse",
+ as.integer(T),
+ as.integer(d),
+ as.integer(as.vector(RVM$family)),
+ as.integer(as.vector(RVM$MaxMat)),
+ as.integer(as.vector(RVM$Matrix)),
+ as.integer(as.vector(RVM$CondDistr$direct)),
+ as.integer(as.vector(RVM$CondDistr$indirect)),
+ as.double(as.vector(RVM$par)),
+ as.double(as.vector(RVM$par2)),
+ as.double(as.vector(data)),
+ as.double(as.vector(hessian)),
+ as.double(as.vector(subhess)),
+ as.double(as.vector(der)),
+ as.double(as.vector(subder)),
+ PACKAGE = 'VineCopula')
+
+ hessian=matrix(out[[11]],dd+tt,dd+tt)
+ subhess=matrix(out[[12]],dd+tt,dd+tt)
+ der=matrix(out[[13]],dd+tt,dd+tt)
+ subder=matrix(out[[14]],dd+tt,dd+tt)
+
+ # der[1:dd,1:dd]=der[dd:1,dd:1]
+ # if(tt>0)
+ # {
+ # der[(dd+1):(dd+tt),1:dd]=der[(dd+tt):(dd+1),dd:1]
+ # der[1:dd,(dd+1):(dd+tt)]=der[dd:1,(dd+tt):(dd+1)]
+ # der[(dd+1):(dd+tt),(dd+1):(dd+tt)]=der[(dd+tt):(dd+1),(dd+tt):(dd+1)]
+ # }
+ # hessian[1:dd,1:dd]=hessian[dd:1,dd:1]
+ # if(tt>0)
+ # {
+ # hessian[(dd+1):(dd+tt),1:dd]=hessian[(dd+tt):(dd+1),dd:1]
+ # hessian[1:dd,(dd+1):(dd+tt)]=hessian[dd:1,(dd+tt):(dd+1)]
+ # hessian[(dd+1):(dd+tt),(dd+1):(dd+tt)]=hessian[(dd+tt):(dd+1),(dd+tt):(dd+1)]
+ # }
+
+ test=apply(hessian,2,function(x) max(abs(x)))
+ hessian=hessian[test>0,test>0]
+ subhess=subhess[test>0,test>0]
+ der=der[test>0,test>0]
+ subder=subder[test>0,test>0]
+
+
+ out=list(hessian=hessian, der=der)
+
+return(out)
+}
Property changes on: pkg/R/RVineHessian.r
___________________________________________________________________
Added: svn:eol-style
+ LF
Modified: pkg/R/RVineLogLik.r
===================================================================
--- pkg/R/RVineLogLik.r 2013-04-04 09:24:02 UTC (rev 5)
+++ pkg/R/RVineLogLik.r 2013-04-04 17:33:10 UTC (rev 6)
@@ -11,7 +11,7 @@
n<-d
N<-T
if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
- if(is(RVM) != "RVineMatrix") stop("'RVM' has to be an RVineMatrix object.")
+ if(!is(RVM, "RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
o = diag(RVM$Matrix)
if(any(o != length(o):1))
Property changes on: pkg/R/RVineLogLik.r
___________________________________________________________________
Added: svn:eol-style
+ LF
Modified: pkg/R/RVineMLE.R
===================================================================
--- pkg/R/RVineMLE.R 2013-04-04 09:24:02 UTC (rev 5)
+++ pkg/R/RVineMLE.R 2013-04-04 17:33:10 UTC (rev 6)
@@ -1,361 +1,361 @@
-RVineMLE <- function(data, RVM, start=RVM$par, start2=RVM$par2, maxit=200, max.df=30, max.BB=list(BB1=c(5,6),BB6=c(6,6),BB7=c(5,6),BB8=c(6,1)), grad=FALSE, hessian=FALSE, se=FALSE, ...)
-{
- if(is(RVM) != "RVineMatrix") stop("'RVM' has to be an RVineMatrix object.")
- if(maxit<=0) stop("'maxit' has to be greater than zero.")
-
- 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.")
- if(max.BB$BB1[1] < 0.001) stop("The upper bound for the first parameter of the BB1 copula should be greater than 0.001 (lower bound for estimation).")
- if(max.BB$BB1[2] < 1.001) stop("The upper bound for the second parameter of the BB1 copula should be greater than 1.001 (lower bound for estimation).")
- if(max.BB$BB6[1] < 1.001) stop("The upper bound for the first parameter of the BB6 copula should be greater than 1.001 (lower bound for estimation).")
- if(max.BB$BB6[2] < 1.001) stop("The upper bound for the second parameter of the BB6 copula should be greater than 1.001 (lower bound for estimation).")
- if(max.BB$BB7[1] < 1.001) stop("The upper bound for the first parameter of the BB7 copula should be greater than 1.001 (lower bound for estimation).")
- if(max.BB$BB7[2] < 0.001) stop("The upper bound for the second parameter of the BB7 copula should be greater than 0.001 (lower bound for estimation).")
- if(max.BB$BB8[1] < 1.001) stop("The upper bound for the first parameter of the BB1 copula should be greater than 0.001 (lower bound for estimation).")
- if(max.BB$BB8[2] < 0.001 || max.BB$BB8[2] > 1) stop("The upper bound for the second parameter of the BB1 copula should be in the interval [0,1].")
-
-
- Matrix=RVM$Matrix
-
- if(!all(start %in% c(0,NA)))
- {
- for(i in 2:dim(Matrix)[1]){
- for(j in 1:(i-1)){
- if((RVM$family[i,j]==1 || RVM$family[i,j]==2) && abs(start[i,j])>=1) stop("The parameter of the Gaussian and t-copula has to be in the interval (-1,1).")
- if(RVM$family[i,j]==2 && start2[i,j]<=1) stop("The degrees of freedom parameter of the t-copula has to be larger than 1.")
- if((RVM$family[i,j]==3 || RVM$family[i,j]==13) && start[i,j]<=0) stop("The parameter of the Clayton copula has to be positive.")
- if((RVM$family[i,j]==4 || RVM$family[i,j]==14) && start[i,j]<1) stop("The parameter of the Gumbel copula has to be in the interval [1,oo).")
- if((RVM$family[i,j]==6 || RVM$family[i,j]==16) && start[i,j]<=1) stop("The parameter of the Joe copula has to be in the interval (1,oo).")
- if(RVM$family[i,j]==5 && start[i,j]==0) stop("The parameter of the Frank copula has to be unequal to 0.")
- if((RVM$family[i,j]==7 || RVM$family[i,j]==17) && start[i,j]<=0) stop("The first parameter of the BB1 copula has to be positive.")
- if((RVM$family[i,j]==7 || RVM$family[i,j]==17) && start2[i,j]<1) stop("The second parameter of the BB1 copula has to be in the interval [1,oo).")
- if((RVM$family[i,j]==8 || RVM$family[i,j]==18) && start[i,j]<=0) stop("The first parameter of the BB6 copula has to be in the interval [1,oo).")
- if((RVM$family[i,j]==8 || RVM$family[i,j]==18) && start2[i,j]<1) stop("The second parameter of the BB6 copula has to be in the interval [1,oo).")
- if((RVM$family[i,j]==9 || RVM$family[i,j]==19) && start[i,j]<1) stop("The first parameter of the BB7 copula has to be in the interval [1,oo).")
- if((RVM$family[i,j]==9 || RVM$family[i,j]==19) && start2[i,j]<=0) stop("The second parameter of the BB7 copula has to be positive.")
- if((RVM$family[i,j]==10 || RVM$family[i,j]==20) && start[i,j]<1) stop("The first parameter of the BB8 copula has to be in the interval [1,oo).")
- if((RVM$family[i,j]==10 || RVM$family[i,j]==20) && (start2[i,j]<=0 || start2[i,j]>1)) stop("The second parameter of the BB8 copula has to be in the interval (0,1].")
- if((RVM$family[i,j]==23 || RVM$family[i,j]==33) && start[i,j]>=0) stop("The parameter of the rotated Clayton copula has to be negative.")
- if((RVM$family[i,j]==24 || RVM$family[i,j]==34) && start[i,j]>-1) stop("The parameter of the rotated Gumbel copula has to be in the interval (-oo,-1].")
- if((RVM$family[i,j]==26 || RVM$family[i,j]==36) && start[i,j]>=-1) stop("The parameter of the rotated Joe copula has to be in the interval (-oo,-1).")
- if((RVM$family[i,j]==27 || RVM$family[i,j]==37) && start[i,j]>=0) stop("The first parameter of the rotated BB1 copula has to be negative.")
- if((RVM$family[i,j]==27 || RVM$family[i,j]==37) && start2[i,j]>-1) stop("The second parameter of the rotated BB1 copula has to be in the interval (-oo,-1].")
- if((RVM$family[i,j]==28 || RVM$family[i,j]==38) && start[i,j]>=0) stop("The first parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
- if((RVM$family[i,j]==28 || RVM$family[i,j]==38) && start2[i,j]>-1) stop("The second parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
- if((RVM$family[i,j]==29 || RVM$family[i,j]==39) && start[i,j]>-1) stop("The first parameter of the rotated BB7 copula has to be in the interval (-oo,-1].")
- if((RVM$family[i,j]==29 || RVM$family[i,j]==39) && start2[i,j]>=0) stop("The second parameter of the rotated BB7 copula has to be negative.")
- if((RVM$family[i,j]==30 || RVM$family[i,j]==40) && start[i,j]>-1) stop("The first parameter of the rotated BB8 copula has to be in the interval (-oo,-1].")
- if((RVM$family[i,j]==30 || RVM$family[i,j]==40) && (start2[i,j]>=0 || start2[i,j]<(-1))) stop("The second parameter of the rotated BB8 copula has to be in the interval [-1,0).")
- }
- }
- }
-
- data<-as.matrix(data)
- if(any(data>1) || any(data<0)) stop("Data has be in the interval [0,1].")
- d=dim(data)[2]
- T=dim(data)[1]
- n<-d
- N<-T
- if(n != dim(RVM)) stop("Dimensions of 'data' and 'RVM' do not match.")
- if(T<2) stop("Number of observations has to be at least 2.")
-
- o = diag(RVM$Matrix)
- oldRVM = RVM
- RVM = normalizeRVineMatrix(RVM)
- data = data[,o[length(o):1]]
-
-
- n = dim(RVM)
- N = dim(data)[1]
-
- if(all(start==0))
- {
- est_start=RVineSeqEst(data,RVM,max.df=max.df,max.BB=max.BB)
- start=est_start$RVM$par
- start2=est_start$RVM$par2
- }
-
- posParams = (RVM$family > 0)
- posParams2 = (RVM$family %in% c(2,7,8,9,10,17,18,19,20,27,28,29,30,37,38,39,40))
-
- posParams[is.na(posParams)] = FALSE
- posParams2[is.na(posParams2)] = FALSE
-
- nParams = sum(posParams, na.rm=TRUE)
- nParams2 = sum(posParams2, na.rm=TRUE)
-
- startpar = double(nParams+nParams2)
-
- Copula.Types = RVM$family[posParams]
-
- startpar[1:nParams] = start[posParams]
- if(nParams2 > 0){
- startpar[(nParams+1):(nParams+nParams2)] = start2[posParams2]
- }
-
- #Grenzen
- lb = double(nParams+nParams2)
- ub = double(nParams+nParams2)
-
- for(i in 1:nParams){
- if(Copula.Types[i] %in% c(1,2)){ #Normal
- lb[i] = -0.98
- ub[i] = 0.98
- }else if(Copula.Types[i] %in% c(3,13)){ #clayton
- lb[i] = 0.0001
- ub[i] = Inf
- }else if(Copula.Types[i] %in% c(23,33)){ #rotated clayton
- lb[i] = -Inf
- ub[i] = -0.0001
- }else if(Copula.Types[i] %in% c(4,14)){ #gumbel
- lb[i] = 1.0001
- ub[i] = Inf
- }else if(Copula.Types[i] %in% c(24,34)){ #rotated gumbel
- lb[i] = -Inf
- ub[i] = -1.0001
- }else if(Copula.Types[i]==5){ #frank
- lb[i] = -Inf
- ub[i] = Inf
- }else if(Copula.Types[i] %in% c(6,16)){ #joe
- lb[i] = 1.0001
- ub[i] = Inf
- }else if(Copula.Types[i] %in% c(26,36)){ #rotated joe
- lb[i] = -Inf
- ub[i] = -1.0001
- }else if(Copula.Types[i] %in% c(7,17)){ #bb1
- lb[i] = 0.001
- ub[i] = max.BB$BB1[1]
- }else if(Copula.Types[i] %in% c(8,18)){ #bb6
- lb[i] = 1.001
- ub[i] = max.BB$BB6[1]
- }else if(Copula.Types[i] %in% c(9,19)){ #bb7
- lb[i] = 1.001
- ub[i] = max.BB$BB7[1]
- }else if(Copula.Types[i] %in% c(10,20)){ #bb8
- lb[i] = 1.001
- ub[i] = max.BB$BB8[1]
- }else if(Copula.Types[i] %in% c(27,37)){ #rotated bb1
- lb[i] = -max.BB$BB1[1]
- ub[i] = -0.001
- }else if(Copula.Types[i] %in% c(28,38)){ #rotated bb6
- lb[i] = -max.BB$BB6[1]
- ub[i] = -1.001
- }else if(Copula.Types[i] %in% c(29,39)){ #rotated bb7
- lb[i] = -max.BB$BB7[1]
- ub[i] = -1.001
- }else if(Copula.Types[i] %in% c(30,40)){ #rotated bb8
- lb[i] = -max.BB$BB8[1]
- ub[i] = -1.001
- }else if(Copula.Types[i] %in% c(43,44)){ #double Clayton, Gumbel
- lb[i] = -0.9999
- ub[i] = 0.9999
- }
-
- }
-
- if(nParams2 > 0){
- todo = which(Copula.Types %in% c(2,7:10,17:20,27:30,37:40))
-
- for(i in 1:nParams2){
- if(Copula.Types[todo[i]]==2){ #t
- lb[nParams+i] = 2.001
- ub[nParams+i] = max.df
- }else if(Copula.Types[todo[i]] %in% c(7,17)){ #bb1
- lb[nParams+i] = 1.001
- ub[nParams+i] = max.BB$BB1[2]
- }else if(Copula.Types[todo[i]] %in% c(8,18)){ #bb6
- lb[nParams+i] = 1.001
- ub[nParams+i] = max.BB$BB6[2]
- }else if(Copula.Types[todo[i]] %in% c(9,19)){ #bb7
- lb[nParams+i] = 0.001
- ub[nParams+i] = max.BB$BB7[2]
- }else if(Copula.Types[todo[i]] %in% c(10,20)){ #bb8
- lb[nParams+i] = 0.001
- ub[nParams+i] = max.BB$BB1[2]
- }else if(Copula.Types[todo[i]] %in% c(27,37)){ #rotated bb1
- lb[nParams+i] = -max.BB$BB1[2]
- ub[nParams+i] = -1.001
- }else if(Copula.Types[todo[i]] %in% c(28,38)){ #rotated bb6
- lb[nParams+i] = -max.BB$BB6[2]
- ub[nParams+i] = -1.001
- }else if(Copula.Types[todo[i]] %in% c(29,39)){ #rotated bb7
- lb[nParams+i] = -max.BB$BB7[2]
- ub[nParams+i] = -1.001
- }else if(Copula.Types[todo[i]] %in% c(30,40)){ #rotated bb8
- lb[nParams+i] = -max.BB$BB8[2]
- ub[nParams+i] = -0.001
- }
- }
- }
-
- startpar1=startpar[Copula.Types!=0]
-
- #calcupdate=NA
- #V=NA
-
- optim_LL <- function(parm,data,posParams,posParams2,Copula.Types,start,start2,RVM, calcupdate=NA){
-
- nParams = sum(posParams, na.rm=TRUE)
- nParams2 = sum(posParams2, na.rm=TRUE)
-
- #for(k in 1:nParams)
- #{
- # if(Copula.Types[k]==0)
- # {
- # if(k==1){parm=c(0,parm)}
- # else if(k>length(parm)){parm=c(parm,0)}
- # else {parm=c(parm[1:(k-1)],0,parm[k:length(parm)])}
- # }
- #}
-
- matrixParams = start
- matrixParams2 = start2
-
- matrixParams[posParams] <- parm[1:nParams]
-
- if(nParams2 > 0){
- matrixParams2[posParams2] <- parm[(nParams+1):(nParams+nParams2)]
- }
-
- ll=RVineLogLik(data,RVM,par=matrixParams, par2=matrixParams2)
-
- #V=ll$V
-
- if(is.finite(ll$loglik))
- {
- return(ll$loglik)
- }
- else
- {
- if(is.na(ll$loglik)){message(parm)}
- message(ll$loglik)
- return(-10^305)
- }
- }
-
-
- ableitung <- function(parm,data,posParams,posParams2,Copula.Types,start,start2,RVM,calcupdate)
- {
- nParams = sum(posParams, na.rm=TRUE)
- nParams2 = sum(posParams2, na.rm=TRUE)
-
- # outparm=parm
- # for (i in 1:length(parm)) {
- # handle_parm=parm
- # handle_parm[i]=handle_parm[i]+0.000001
- # handle_parm2=parm
- # handle_parm2[i]=handle_parm2[i]-0.000001
- # outparm[i]=(optim_LL(handle_parm,data,posParams,posParams2,Copula.Types,start,start2,RVM,calcupdate=NA)-optim_LL(handle_parm2,data,posParams,posParams2,Copula.Types,start,start2,RVM,calcupdate=NA))/(2*0.000001)
- # }
- # print("finite differences:")
- # print(outparm)
-
- #for(k in 1:nParams)
- #{
- # if(Copula.Types[k]==0)
- # {
- # if(k==1){parm=c(0,parm)}
- # else if(k>length(parm)){parm=c(parm,0)}
- # else {parm=c(parm[1:(k-1)],0,parm[k:length(parm)])}
- # }
- # }
-
- matrixParams = start
- matrixParams2 = start2
-
- matrixParams[posParams] <- parm[1:nParams]
- if(nParams2 > 0){
- matrixParams2[posParams2] <- parm[(nParams+1):(nParams+nParams2)]
- }
-
- grad=RVineGrad(data=data,RVM=RVM,par=matrixParams, par2=matrixParams2, posParams=posParams)$gradient
- return(grad)
- }
-
- pscale=numeric()
- for(i in 1:nParams)
- {
- pscale[i] = ifelse(Copula.Types[i] %in% c(1,2,43,44), 0.01, 1)
- }
- pscale=c(pscale,rep(1,nParams2))
-
- if(!exists("factr")) # Toleranz etwas hoch setzen (gröber)
- factr=1e8
-
- if(all(Copula.Types %in% c(0,1,2,3:6,13,14,16,23,24,26,33,34,36,43,44)) && grad==TRUE)
- {
-
- #n=dim(RVM)
- #calcupdate=array(0,dim=c(n,n,n,n))
- #for(i in (n-1):1){
- # for(k in n:(i+1)){
- # calcupdate[, ,k,i ]=RVineMatrixUpdate(RVM,k,i)
- # }
- #}
- if(hessian==TRUE || se==TRUE)
- {
- out1 = optim(par=startpar1,fn=optim_LL, gr=ableitung,
- data=data,posParams=posParams,posParams2=posParams2,Copula.Types=Copula.Types,start=start,start2=start2,RVM=RVM,
- method="L-BFGS-B",control=list(fnscale=-1,maxit=maxit,trace=1,parscale=pscale,factr=factr),lower=lb,upper=ub,hessian=TRUE)
- }
- else
- {
- out1 = optim(par=startpar1,fn=optim_LL, gr=ableitung,
- data=data,posParams=posParams,posParams2=posParams2,Copula.Types=Copula.Types,start=start,start2=start2,RVM=RVM,
- method="L-BFGS-B",control=list(fnscale=-1,maxit=maxit,trace=1,parscale=pscale,factr=factr),lower=lb,upper=ub)
- }
- }
- else
- {
- if(hessian==TRUE || se==TRUE)
- {
- out1 = optim(par=startpar1,fn=optim_LL,
- data=data,posParams=posParams,posParams2=posParams2,Copula.Types=Copula.Types,start=start,start2=start2,RVM=RVM, calcupdate=NA,
- method="L-BFGS-B",control=list(fnscale=-1,maxit=maxit,trace=1,parscale=pscale,factr=factr,...),lower=lb,upper=ub,hessian=TRUE)
- }
- else
- {
- out1 = optim(par=startpar1,fn=optim_LL,
- data=data,posParams=posParams,posParams2=posParams2,Copula.Types=Copula.Types,start=start,start2=start2,RVM=RVM, calcupdate=NA,
- method="L-BFGS-B",control=list(fnscale=-1,maxit=maxit,trace=1,parscale=pscale,factr=factr,...),lower=lb,upper=ub)
- }
- }
-
- out = list()
-
- out$value = out1$value
- out$convergence = out1$convergence
- out$message = out1$message
- out$counts=out1$counts
- if(hessian==TRUE)
- out$hessian = out1$hessian
-
- if(se==TRUE)
- out1$se=sqrt((diag(solve(-out1$hessian))))
-
- out$RVM = oldRVM
-
- kk=1
- for (ll in 1:nParams)
- {
- out1$par[ll]=out1$par[ll]
- if(Copula.Types[ll] %in% c(2,7:10,17:20,27:30,37:40))
- {
- out1$par[nParams+kk]=out1$par[nParams+kk]
- kk=kk+1
- }
- }
-
- out$RVM$par[posParams] = out1$par[1:nParams]
- out$RVM$par2[posParams2] = out1$par[(nParams+1):(nParams+nParams2)]
-
- if(se==TRUE)
- {
- out$se=matrix(0,d,d)
- out$se2=matrix(0,d,d)
- out$se[posParams] = out1$se[1:nParams]
- out$se2[posParams2] = out1$se[(nParams+1):(nParams+nParams2)]
- }
-
- return(out)
-}
+RVineMLE <- function(data, RVM, start=RVM$par, start2=RVM$par2, maxit=200, max.df=30, max.BB=list(BB1=c(5,6),BB6=c(6,6),BB7=c(5,6),BB8=c(6,1)), grad=FALSE, hessian=FALSE, se=FALSE, ...)
+{
+ if(!is(RVM, "RVineMatrix")) stop("'RVM' has to be an RVineMatrix object.")
+ if(maxit<=0) stop("'maxit' has to be greater than zero.")
+
+ 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.")
+ if(max.BB$BB1[1] < 0.001) stop("The upper bound for the first parameter of the BB1 copula should be greater than 0.001 (lower bound for estimation).")
+ if(max.BB$BB1[2] < 1.001) stop("The upper bound for the second parameter of the BB1 copula should be greater than 1.001 (lower bound for estimation).")
+ if(max.BB$BB6[1] < 1.001) stop("The upper bound for the first parameter of the BB6 copula should be greater than 1.001 (lower bound for estimation).")
+ if(max.BB$BB6[2] < 1.001) stop("The upper bound for the second parameter of the BB6 copula should be greater than 1.001 (lower bound for estimation).")
+ if(max.BB$BB7[1] < 1.001) stop("The upper bound for the first parameter of the BB7 copula should be greater than 1.001 (lower bound for estimation).")
+ if(max.BB$BB7[2] < 0.001) stop("The upper bound for the second parameter of the BB7 copula should be greater than 0.001 (lower bound for estimation).")
+ if(max.BB$BB8[1] < 1.001) stop("The upper bound for the first parameter of the BB1 copula should be greater than 0.001 (lower bound for estimation).")
+ if(max.BB$BB8[2] < 0.001 || max.BB$BB8[2] > 1) stop("The upper bound for the second parameter of the BB1 copula should be in the interval [0,1].")
+
+
+ Matrix=RVM$Matrix
+
+ if(!all(start %in% c(0,NA)))
+ {
+ for(i in 2:dim(Matrix)[1]){
+ for(j in 1:(i-1)){
+ if((RVM$family[i,j]==1 || RVM$family[i,j]==2) && abs(start[i,j])>=1) stop("The parameter of the Gaussian and t-copula has to be in the interval (-1,1).")
+ if(RVM$family[i,j]==2 && start2[i,j]<=1) stop("The degrees of freedom parameter of the t-copula has to be larger than 1.")
+ if((RVM$family[i,j]==3 || RVM$family[i,j]==13) && start[i,j]<=0) stop("The parameter of the Clayton copula has to be positive.")
+ if((RVM$family[i,j]==4 || RVM$family[i,j]==14) && start[i,j]<1) stop("The parameter of the Gumbel copula has to be in the interval [1,oo).")
+ if((RVM$family[i,j]==6 || RVM$family[i,j]==16) && start[i,j]<=1) stop("The parameter of the Joe copula has to be in the interval (1,oo).")
+ if(RVM$family[i,j]==5 && start[i,j]==0) stop("The parameter of the Frank copula has to be unequal to 0.")
+ if((RVM$family[i,j]==7 || RVM$family[i,j]==17) && start[i,j]<=0) stop("The first parameter of the BB1 copula has to be positive.")
+ if((RVM$family[i,j]==7 || RVM$family[i,j]==17) && start2[i,j]<1) stop("The second parameter of the BB1 copula has to be in the interval [1,oo).")
+ if((RVM$family[i,j]==8 || RVM$family[i,j]==18) && start[i,j]<=0) stop("The first parameter of the BB6 copula has to be in the interval [1,oo).")
+ if((RVM$family[i,j]==8 || RVM$family[i,j]==18) && start2[i,j]<1) stop("The second parameter of the BB6 copula has to be in the interval [1,oo).")
+ if((RVM$family[i,j]==9 || RVM$family[i,j]==19) && start[i,j]<1) stop("The first parameter of the BB7 copula has to be in the interval [1,oo).")
+ if((RVM$family[i,j]==9 || RVM$family[i,j]==19) && start2[i,j]<=0) stop("The second parameter of the BB7 copula has to be positive.")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/vinecopula -r 6
Mehr Informationen über die Mailingliste Vinecopula-commits