Laura Gianfagna
2015-Mar-31 16:39 UTC
[R] Does fitCopula work for amhCopula and joeCopula?
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]]
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.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]] > > ______________________________________________ > R-help at 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.mes at cbs.dk Priv: PDalgd at gmail.com