Does fitCopula work for amhCopula and joeCopula?

On 31 Mar 2015, at 18:39 , Laura Gianfagna <laura.gianfagna at outlook.it> wrote:

Good evening, this is a part of my Routine  which calculates the copula parameter and loglikelihood for each pair of rows of a data matrix, choosing, for each pair, the copula which gives the maximum likelihood. If I do my computation with this routine with only:


f <- frankCopula(2,2)
 g <- gumbelCopula(2,2)
 c <- claytonCopula(2,2)


the program works correctly and gives the expected results.

If  I insert also:


 a <- amhCopula(1,2)
 j <- joeCopula(2,2)


then the program doesn?t work anymore. 

I tried on samples such as:


n <- 1000
f <- frankCopula(20,2)
x_1 <- rCopula(n,f)
f <- gumbelCopula(50,2)
x_2 <- rCopula(n,f)
f <- joeCopula(70,2)
x_3<- rCopula(n,f)
x <- cbind(x_1, x_2, x_3)
data <- t(x)
dim <- dim(data)[1]





Here is the part of code of Routine_Copula:

Routine_Copula <- function(data,dim){

 library(copula)
 library(gtools)

 n <- dim(data)[1];  # number of rows of the input matrix
 m <- dim(data)[2];  # number of columns of the input matrix

 # Probability integral transform of the data
 ecdf <- matrix(0,n,m);
 for (i in 1:n){
   e <- matrix(data[i,],m,1);
   #ecdf[i,] <- pobs(e);
   ecdf[i,] <- pobs(e, na.last=TRUE);
   #na.last for controlling the treatment of NAs. If TRUE, missing values in the data are put last; if FALSE, they are put first; if NA, they are removed; if "keep" they are kept with rank NA.

 }



f <- frankCopula(2,2)
 g <- gumbelCopula(2,2)
 c <- claytonCopula(2,2)
 a <- amhCopula(1,2)
 j <- joeCopula(2,2)





[?.]


for (j in 1:n_comb){
   input <- t(ecdf[comb[,j],])

   try(summary <- fitCopula(f,input,method='mpl',start=2),silent=TRUE);
   resmatpar[j,1] <- summary at estimate;
   resmatllk[j,1] <- summary at loglik;

   try(summary <- fitCopula(g,input,method='mpl',start=2),silent=TRUE);
   resmatpar[j,2] <- summary at estimate;
   resmatllk[j,2] <- summary at loglik;

   try(summary <- fitCopula(c,input,method='mpl',start=2),silent=TRUE);
   resmatpar[j,3] <- summary at estimate;
   resmatllk[j,3] <- summary at loglik;


try(summary <- fitCopula(a,input,method='mpl',start=1),silent=TRUE);
   resmatpar[j,4] <- summary at estimate;
    resmatllk[j,4] <- summary at loglik;

try(summary <- fitCopula(j,input,method='mpl',start=2),silent=TRUE);     

resmatpar[j,5] <- summary at estimate;
resmatllk[j,5] <- summary at loglik;

d <- c(resmatllk[j,1],resmatllk[j,2],resmatllk[j,3],resmatllk[j,4],resmatllk[j,5]);



   copchoice[j] <- which(d==max(d));
   param[j] <- resmatpar[j,copchoice[j]];
   loglik[j] <- resmatllk[j,copchoice[j]];

 }


Thank you

Laura Gianfagna


?
	[[alternative HTML version deleted]]