You might be better off on a list specifically for financial data analysis. In
particular, few people here will instantly know that you are talking about
functions the copula package. And even many of us who have dabbled in copulas
will have difficulty knowing all of the details of them. Also, "doesn't work"
is too unspecific to allow people to pinpoint the problem, and you might be
expected to narrow down and debug the issue rather than just handing us a lot
of R code.
Offhand, you certainly have a problem with reuse of the variable name `j`.
-pd
> On 31 Mar 2015, at 18:39 , Laura Gianfagna <laura.gianfa...@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@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.
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd....@cbs.dk Priv: pda...@gmail.com
______________________________________________
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.
