[R] Does fitCopula work for amhCopula and joeCopula?
Laura Gianfagna
laura.gianfagna at outlook.it
Tue Mar 31 18:39:59 CEST 2015
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]]
More information about the R-help
mailing list