[R] Does fitCopula work for amhCopula and joeCopula?

[Laura Gianfagna]

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@estimate;
    resmatllk[j,1] <- summary@loglik;
    
    try(summary <- fitCopula(g,input,method='mpl',start=2),silent=TRUE);
    resmatpar[j,2] <- summary@estimate;
    resmatllk[j,2] <- summary@loglik;
    
    try(summary <- fitCopula(c,input,method='mpl',start=2),silent=TRUE);
    resmatpar[j,3] <- summary@estimate;
    resmatllk[j,3] <- summary@loglik;


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

 resmatpar[j,5] <- summary@estimate;
resmatllk[j,5] <- summary@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]]

______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.