R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## Copyright (C) 2021 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan > ## > ## This program is free software; you can redistribute it and/or modify it under > ## the terms of the GNU General Public License as published by the Free Software > ## Foundation; either version 3 of the License, or (at your option) any later > ## version. > ## > ## This program is distributed in the hope that it will be useful, but WITHOUT > ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS > ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more > ## details. > ## > ## You should have received a copy of the GNU General Public License along with > ## this program; if not, see . > > > require(copula) Loading required package: copula > source(system.file("Rsource", "tstFit-fn.R", package="copula", mustWork=TRUE)) > source(system.file("Rsource", "utils.R", package="copula", mustWork=TRUE)) Loading required package: tools > ##-> assertError(), assert.EQ(), ... showProc.time() > > (doExtras <- copula:::doExtras()) [1] FALSE > > > uu <- array(c(9, 7, 8, 3, 2, 4, 1, 5, 6, 10, + 6, 9, 1, 7, 3, 2, 5, 8, 4, 10), dim = c(10L, 2L)) / 11 > suppressWarnings(RNGversion("3.5.0")) ## << FIXME: eliminate - adapting where needed below ! > set.seed(7) > u3 <- cbind(uu, round(runif(10),2)) > > ### t-copula instead of normal -- minimal set for testing here: > ## d = 2 > ## fit1() here catches the error: "irho" not available for "tCopula": > (f1 <- fit1(tCopula(df.fixed=TRUE), x = uu)) Caught fitCopula() error: 'unable to find an inherited method for function 'iRho' for signature 'copula = "tCopula"'' *-*-* est se itau 0.3746066 0.3867052 irho NA NA mpl 0.4551169 0.4190495 ml 0.4551169 0.3094090 > stopifnot(identical(f1, fit1(tCopula(df.fixed=TRUE), + x = data.frame(uu))))# *WITH* a warning; keep data.frame() working! Caught fitCopula() error: 'unable to find an inherited method for function 'iRho' for signature 'copula = "tCopula"'' *-*-* Warning message: In .local(copula, data, ...) : coercing 'data' to a matrix. > ## did not work with data.frame before 2012-08-12 > > ## for df.fixed=FALSE, have 2 parameters ==> cannot use "fit1": > (f2.t <- fitCopula(tCopula(), uu, method="itau"))# (+ warning: ... 'df.fixed=TRUE' ) Call: fitCopula(tCopula(), data = uu, ... = pairlist(method = "itau")) Fit based on "inversion of Kendall's tau" and 10 2-dimensional observations. Copula: tCopula rho.1 0.3746 Warning message: In fitCopula.icor(copula, x = data, method = method, estimate.variance = estimate.variance, : "itau" fitting ==> copula coerced to 'df.fixed=TRUE' > assertError( ## 14.Jan.2016: irho() has been "non-sense" -> now erronous: + fitCopula(tCopula(), uu, method="irho"), verbose=TRUE) Asserted error: unable to find an inherited method for function 'iRho' for signature 'copula = "tCopula"' > showProc.time() Time (user system elapsed): 0.2 0.02 0.22 > if(doExtras) { + print(f2.m <- fitCopula(tCopula(), uu, method= "ml")) + print(f2.M <- fitCopula(tCopula(), uu, method= "mpl")) + print(summary(f2.m)) # gives SE for 'df' {from optim()} + print(summary(f2.M)) # no SE for 'df' (for now ..) + stopifnot(all.equal(coef(f2.m), coef(f2.M))) + showProc.time() + } > > ## d = 3 : ------------- > ## ok with df.fixed > tC3f <- tCopula(dim=3, dispstr="un", df.fixed=TRUE) > tC3 <- tCopula(dim=3, dispstr="un") > f3 <- fitCopula(tC3f, u3, method="itau") > f3.t <- fitCopula(tC3 , u3, method="itau") # warning: coercing to df.fixed=TRUE Warning message: In fitCopula.icor(copula, x = data, method = method, estimate.variance = estimate.variance, : "itau" fitting ==> copula coerced to 'df.fixed=TRUE' > summary(f3.t) Call: fitCopula(tC3, data = u3, ... = pairlist(method = "itau")) Fit based on "inversion of Kendall's tau" and 10 3-dimensional observations. t-copula, dim. d = 3 Estimate Std. Error rho.1 0.3746 0.387 rho.2 0.3090 0.326 rho.3 0.3746 0.405 > cf3 <- coef(f3, SE = TRUE) > stopifnot(all.equal(unname(cf3), cbind(c(0.374607, 0.309017, 0.374607), + c(0.386705, 0.325995, 0.405493)), + tol = 5e-5), # seen 6e-7 + all.equal(coef(f3), coef(f3.t))) > if(FALSE) ## Error: iRho() method for class "tCopula" not yet implemented + (f3.r <- fitCopula(tC3, u3, method="irho")) > showProc.time() Time (user system elapsed): 0.02 0 0.01 > > octr <- list(trace=TRUE, REPORT=1, reltol = 1e-5)# less accurate than by default > if(doExtras) { + fitC.m <- function(C, data, method, ...) + fitCopula(C, data, method=method, optim.control=octr, traceOpt=TRUE, ...) + print(f3.m <- fitC.m(tC3, u3, method= "ml")); c.m <- coef(f3.m, SE=TRUE) + print(f3.M <- fitC.m(tC3, u3, method= "mpl")); c.M <- coef(f3.M, SE=TRUE) + showProc.time() + is.matrix(print(c.m)) + is.character(SE <- "Std. Error") + stopifnot(all.equal(c.m[,1], c.M[,1]), # the estimates don't differ; the SE's do: + c.M[-4,SE] >= c.m[-4,SE], + c.m[, SE] > 0) + } > > set.seed(17) > d <- 5 # dimension > nu <- 4 # degrees of freedom > ## define and sample the copula, build pseudo-observations > ec4 <- tCopula(dim=d, df=nu, df.fixed=TRUE) # <- copula with param NA > (r <- iTau(ec4, tau <- c(0.2, 0.4, 0.6))) [1] 0.3090170 0.5877853 0.8090170 > P <- c(r[2], r[1], r[1], r[1], # upper triangle (w/o diagonal) of corr.matrix + r[1], r[1], r[1], + r[3], r[3], + r[3]) > assertError( setTheta(ec4, value = P) ) > validObject(ex4. <- setTheta(ec4, value = 0.8)) # set the only (non-fixed) parameter [1] TRUE > ## TODO "check" getTheta(ex4., ...) > ## rather need "un" dispersion: Now with smarter tCopula(): > (uc4 <- tCopula(dim=d, df=nu, disp = "un", df.fixed=FALSE)) t-copula, dim. d = 5 Dimension: 5 Parameters: rho.1 = NA rho.2 = NA rho.3 = NA rho.4 = NA rho.5 = NA rho.6 = NA rho.7 = NA rho.8 = NA rho.9 = NA rho.10 = NA df = 4 dispstr: un > validObject(uc4p <- setTheta(uc4, value = c(P, df=nu))) [1] TRUE > U. <- pobs(rCopula(n=1000, copula=uc4p)) > splom2(U.) # => now correct dependency > (cU <- cor(U., method="kendall")) # => correct: [,1] [,2] [,3] [,4] [,5] [1,] 1.0000000 0.4271471 0.2095616 0.2245806 0.1969489 [2,] 0.4271471 1.0000000 0.2114194 0.2064505 0.2156076 [3,] 0.2095616 0.2114194 1.0000000 0.6196877 0.6062623 [4,] 0.2245806 0.2064505 0.6196877 1.0000000 0.6062422 [5,] 0.1969489 0.2156076 0.6062623 0.6062422 1.0000000 > stopifnot(cor(P, cU[lower.tri(cU)]) > 0.99) > > ## Fitting a t-copula with "itau.mpl" with disp="un" > (fm4u <- fitCopula(uc4, U., method="itau.mpl", traceOpt = TRUE)) 1/nu= 12.2229124 => nu= 0.0818; logL= -35081.013 1/nu= 19.7770876 => nu= 0.0506; logL= -62822.427 1/nu= 7.55417528 => nu= 0.1324; logL= -18580.952 1/nu= 4.66873708 => nu= 0.2142; logL= -8976.378 1/nu= 2.8854382 => nu= 0.3466; logL= -3558.4793 1/nu= 1.78329888 => nu= 0.5608; logL= -640.42313 1/nu= 1.10213932 => nu= 0.9073; logL= 819.01527 1/nu= 0.68115956 => nu= 1.4681; logL= 1462.6807 1/nu= 0.42097976 => nu= 2.3754; logL= 1689.1506 1/nu= 0.100898586 => nu= 9.9109; logL= 1677.214 1/nu= 0.29871963 => nu= 3.3476; logL= 1730.6775 1/nu= 0.270724479 => nu= 3.6938; logL= 1733.1412 1/nu= 0.258450766 => nu= 3.8692; logL= 1733.3212 1/nu= 0.260560314 => nu= 3.8379; logL= 1733.3302 1/nu= 0.260650866 => nu= 3.8365; logL= 1733.3302 1/nu= 0.260610172 => nu= 3.8371; logL= 1733.3302 1/nu= 0.26069156 => nu= 3.8360; logL= 1733.3302 1/nu= 0.260650866 => nu= 3.8365; logL= 1733.3302 Call: fitCopula(uc4, data = U., ... = pairlist(method = "itau.mpl", traceOpt = TRUE)) Fit based on "itau for dispersion matrix P and maximum likelihood for df" and 1000 5-dimensional observations. Copula: tCopula rho.1 rho.2 rho.3 rho.4 rho.5 rho.6 rho.7 rho.8 rho.9 rho.10 df 0.6217 0.3233 0.3455 0.3045 0.3260 0.3186 0.3322 0.8268 0.8148 0.8147 3.8365 The maximized loglikelihood is 1733 Optimization converged > ## Fitting t-copulas ............. with disp = "ex" and "ar" : > uc4.ex <- tCopula(dim=d, df=nu, disp = "ex", df.fixed=FALSE) > uc4.ar <- tCopula(dim=d, df=nu, disp = "ar1", df.fixed=FALSE) > validObject(uc4p.ex <- setTheta(uc4.ex, value = c(0.75, df=nu))) [1] TRUE > validObject(uc4p.ar <- setTheta(uc4.ar, value = c(0.75, df=nu))) [1] TRUE > U.ex <- pobs(rCopula(n=1000, copula=uc4p.ex)) > U.ar <- pobs(rCopula(n=1000, copula=uc4p.ar)) > if(FALSE) { # The following are not available (yet); see ~/R/fitCopula.R + ## Fitting a t-copula with "itau.mpl" with disp="ex" + (fm4e <- fitCopula(uc4.ex, U.ex, method="itau.mpl")) + ## Fitting a t-copula with "itau.mpl" with disp="ar" + (fm4e <- fitCopula(uc4.ar, U.ar, method="itau.mpl")) + } > > ## Extra checks -------------------------------------------------------- > if(doExtras) { ## Typically not run on R CMD check + tCop <- tCopula(c(0.2,0.4,0.6), dim=3, dispstr="un", df=5) + set.seed(101) + x <- rCopula(n=200, tCop) # "true" observations (simulated) + ## Maximum likelihood (start = (rho[1:3], df)) + print(summary(tc.ml <- + fitCopula(tCopula(dim=3, dispstr="un"), x, method="ml", + start=c(0,0,0, 10)))) + print(summary(tc.ml. <- + fitCopula(tCopula(dim=3, dispstr="un"), x, method="ml")))# w/o 'start' + + ## Maximum pseudo-likelihood (the asymptotic variance cannot be estimated) + u <- pobs(x) + print(tc.mpl <- fitCopula(tCopula(dim=3, dispstr="un"), + u, method="mpl", estimate.variance=FALSE, + start= c(0,0,0,10))) + ## Without 'start' + tc.mp. <- fitCopula(tCopula(dim=3, dispstr="un"), x, estimate.variance=FALSE) + print(tc.mp.) + noC <- function(x) { x@fitting.stats$counts <- NULL; x@call <- quote(dummy()); x } + + assert.EQ(noC(tc.ml) , noC(tc.ml.), tol= .005) # nothing + assert.EQ(noC(tc.mpl), noC(tc.mp.), tol= .100, giveRE=TRUE) # shows diff + + ## The same t copula but with df.fixed=TRUE (=> use the same data!) + tC3u5 <- tCopula(dim=3, dispstr="un", df=5, df.fixed=TRUE) + ## Maximum likelihood (start = rho[1:3]) + print(tcF.ml <- fitCopula(tC3u5, x, method="ml", start=c(0,0,0))) + print(tcF.ml. <- fitCopula(tC3u5, x, method="ml")) # without 'start' + assert.EQ(noC(tcF.ml), noC(tcF.ml.), tol= 4e-4) + print(vcov(tcF.ml)) # the (estimated, asymptotic) var-cov matrix + + ## Maximum pseudo-likelihood (the asymptotic variance cannot be estimated) + print(tcF.mpl <- fitCopula(tC3u5, u, method="mpl", estimate.variance=FALSE, + start=c(0,0,0))) + print(tcF.mp. <- fitCopula(tC3u5, u, method="mpl", estimate.variance=FALSE)) + assert.EQ(noC(tcF.mpl), noC(tcF.mp.), tol = 1e-5) + } # end Extras ---------------------------------------------------- > > ## fitMvdc() -- first 2 D -- from Yiyun Shou's bug report: --------------------- > > ct.2 <- tCopula(param=0.2, dim=2, dispstr = "ex", df.fixed=TRUE) > mvt.2.ne <- mvdc(copula = ct.2, margins = c("norm", "exp"), + paramMargins = list(list(mean = 0, sd = 2), list(rate = 2))) > mvt.2.ne ## --> four free parameters in total: rho, mean, sd, and rate: Multivariate Distribution Copula based ("mvdc") @ copula: t-copula, dim. d = 2 Dimension: 2 Parameters (partly fixed, see ':='): rho.1 = 0.2 df := 4.0 @ margins: [1] "norm" "exp" with 2 (not identical) margins; with parameters (@ paramMargins) List of 2 $ :List of 2 ..$ mean: num 0 ..$ sd : num 2 $ :List of 1 ..$ rate: num 2 > > getTheta(mvt.2.ne, attr = TRUE) # - shows only the copula parameter -- FIXME ! rho.1 0.2 attr(,"param.lowbnd") [1] -1 attr(,"param.upbnd") [1] 1 > > ## simulate data and fit: > set.seed(17); x.samp <- rMvdc(250, mvt.2.ne) > fit2ne <- fitMvdc(x.samp, mvt.2.ne, start= c(1,1,1, rho = 0.5), + optim.control = list(trace = TRUE), hideWarnings=FALSE) initial value 1126.329900 iter 10 value 600.190021 final value 599.478902 converged initial value 599.478902 final value 599.478902 stopped after 1 iterations Warning messages: 1: In pnorm(x, mean = -15.6081823581344, sd = -1.26949880219363) : NaNs produced 2: In dnorm(x, mean = -15.6081823581344, sd = -1.26949880219363) : NaNs produced 3: In pexp(x, rate = -0.813186200556062) : NaNs produced 4: In dexp(x, rate = -0.813186200556062) : NaNs produced > summary(fit2ne) Call: fitMvdc(data = x.samp, mvdc = mvt.2.ne, start = c(1, 1, 1, rho = 0.5), optim.control = list(trace = TRUE), hideWarnings = FALSE) Maximum Likelihood estimation based on 250 2-dimensional observations. Copula: tCopula Margin 1 : Estimate Std. Error m1.mean 0.06136 0.125 m1.sd 2.07694 0.087 Margin 2 : Estimate Std. Error m2.rate 2.044 0.123 t-copula, dim. d = 2 Estimate Std. Error rho.1 0.1507 0.071 The maximized loglikelihood is -599.5 Optimization converged Number of loglikelihood evaluations: function gradient 54 17 > (confint(fit2ne) -> ci.2ne) 2.5 % 97.5 % m1.mean -0.18308752 0.3058066 m1.sd 1.90644157 2.2474430 m2.rate 1.80186646 2.2857210 rho.1 0.01228649 0.2891935 > stopifnot(exprs = { + all.equal(coef(fit2ne), + c(m1.mean=0.061359521, m1.sd=2.0769423, + m2.rate=2.0437937, rho.1 = 0.15074002), tol=1e-7)# seen 1.48e-8 + all.equal(c(ci.2ne), + c(-0.18309, 1.9064, 1.8019, 0.012286, + 0.30581, 2.2474, 2.2857, 0.28919), tol = 4e-4) # seen 1.65e-5 + }) > > > > ### Fitting multivariate incl margins --- mvdc -------------------------------------- > ### =========================================== > > set.seed(121) > gumbelC <- gumbelCopula(3, dim=2) > gMvGam <- mvdc(gumbelC, c("gamma","gamma"), param = list(list(2,3), list(4,1))) > gMvGam # now nicely show()s -- the *AUTO*-constructed parameter names Multivariate Distribution Copula based ("mvdc") @ copula: Gumbel copula, dim. d = 2 Dimension: 2 Parameters: alpha = 3 @ margins: [1] "gamma" "gamma" with 2 (not identical) margins; with parameters (@ paramMargins) List of 2 $ :List of 2 ..$ shape: num 2 ..$ rate : num 3 $ :List of 2 ..$ shape: num 4 ..$ rate : num 1 > stopifnot(identical(gMvGam@paramMargins, + list(list(shape = 2, rate = 3), + list(shape = 4, rate = 1)))) > X <- rMvdc(16000, gMvGam) > smoothScatter(X, main = "rMvdc(1600, gMvGam)") > > persp (gMvGam, dMvdc, xlim = c(0,4), ylim=c(0,8)) ## almost discrete ???? > contour(gMvGam, dMvdc, xlim = c(0,2), ylim=c(0,8)) > points(X, pch = ".", cex = 2, col=adjustcolor("blue", 0.5)) > > if(doExtras) { + st <- system.time( + fMv <- fitMvdc(X, gMvGam, start = c(1,1,1,1, 1.3),# method="BFGS", + optim.control= list(trace=TRUE))) + print(st) # ~ 59 sec. (lynne 2015); ~ 47 sec (2020) + print(summary(fMv)) + cfM <- coef(fMv, SE=TRUE) + print( all.equal(c(2,3, 4,1, 3), unname(cfM[,"Estimate"]))) # 0.01518 + stopifnot(all.equal(c(2,3, 4,1, 3), unname(cfM[,"Estimate"]), tol = 0.02)) + } > > pFoo <- function(x, lower.tail=TRUE, log.p=FALSE) + pnorm((x - 5)/20, lower.tail=lower.tail, log.p=log.p) > dFoo <- function(x, lower.tail=TRUE, log.p=FALSE) + 1/20* dnorm((x - 5)/20, lower.tail=lower.tail, log.p=log.p) > qFoo <- qunif # must exist; not used for fitting > > ## 'Foo' distribution has *no* parameters: > mv1 <- mvdc(gumbelC, c("gamma","Foo"), param= list(list(3,1), list())) > validObject(mv1) [1] TRUE > stopifnot(nrow(R <- rMvdc(3, mv1)) == 3, ncol(R) == 2) > ## a wrong way: > assertError( + mvW <- mvdc(gumbelC, c("gamma","Foo"), param= list(list(3,1), list(NULL))) + , verbose=TRUE) Asserted error: invalid class "mvdc" object: 'paramMargins[[2]]' must be named properly > ## >> invalid class “mvdc” object: 'paramMargins[[2]]' must be named properly > > > showProc.time() Time (user system elapsed): 2.65 0.03 2.69 > > ## An example which fails (and should)? -- > ## From: Suzanne Li @...queensu.ca, Date: Fri, 9 Aug 2013 > gumbel <- archmCopula(family = "gumbel",dim = 2) > > set.seed(47)# MM {not important, but still want sure reproducibility} > u <- cbind(runif(10),runif(10)) > cor(u[,1], u[,2], method="kendall") [1] -0.02222222 > ## [1] -0.02222222 -- slightly negative > ## this now gives an error: > try(fgu <- fitCopula(gumbel, u, method = "ml")) Warning messages: 1: In iTauGumbelCopula(copula, tau) : For the Gumbel copula, tau must be >= 0. Replacing negative values by 0. 2: In fitCopula.ml(copula, u = data, method = method, start = start, : optim(*, hessian=TRUE) failed: non-finite finite-difference value [1] > ## Error in optim(start, loglikCopula, lower = lower, upper = upper, method = method, : > ## non-finite finite-difference value [1] > ## In addition: Warning message: > ## In .local(copula, tau, ...) : tau is out of the range [0, 1] > copGumbel@paraInterval # -> [1, Inf) = exp([0, Inf)) 'interval' object [1, Inf) > length(par <- 2^c((0:32)/16, 2+(1:10)/8)) # 43 [1] 43 > (t1 <- system.time( + llg <- sapply(par, function(p) loglikCopula(param=p, u=u, copula=gumbel)) + )) user system elapsed 0.03 0.00 0.03 > (t2 <- system.time( + llg. <- loglikCopulaMany(par, u=u, copula=gumbel) + )) user system elapsed 0.03 0.00 0.03 > all.equal(llg, llg., tol=0) [1] TRUE > if(dev.interactive()) plot(par, llg, type="b", col=2) > stopifnot(exprs = { + diff(llg) < 0 # so the maximum is for par = 2^0 = 1 --> at *boundary* of interval + all.equal(llg, llg., tol=4e-16) # see 0 diff (Lnx 64b) + }) > > ## FIXME -- "ml" should return the boundary case, or a much better error message > ## These work (with a warning {which is interesting, but maybe should not be a warning} > ## "perfectly": They give the correct boundary case: > fg.itau <- fitCopula(gumbel, u, method = "itau") Warning message: In iTauGumbelCopula(copula, tau) : For the Gumbel copula, tau must be >= 0. Replacing negative values by 0. > fg.irho <- fitCopula(gumbel, u, method = "irho") Warning message: In iRhoGumbelCopula(copula, rho) : For the Gumbel copula, rho must be >= 0. Replacing negative values by 0. > > ## Is it just this problem? > ## well, the likelihood was not ok for large param=theta; now is: > lrgP <- 100*seq(8, 100, by = 3) > llrg <- vapply(lrgP, function(P) loglikCopula(param=P, u=u, copula=gumbel), NA_real_) > stopifnot(is.finite(llrg), diff(llrg) < 0, llrg < -11990)## no longer NaN > if(FALSE) + plot(lrgP, llrg, type="b", col=2) ## all fine > > ## Now is it because we *really* should use elme() and the "nacopula" families? > ## No, this fails too: "outside interval" error {instead of just -Inf !}: > (copG <- onacopulaL("Gumbel", list(NaN,1:2))) Nested Archimedean copula ("outer_nacopula" of dim. 2), with slot 'comp' = (1, 2) and root 'copula' = Archimedean copula ("acopula"), family "Gumbel" and *no* child copulas > ## Estimation -> error for now -- (FIXME!) > try(efm <- emle(u, copG)) Error in cop@copula@dacopula(u, theta, n.MC = n.MC, log = TRUE) : theta is outside its defined interval In addition: Warning message: In initOpt(cop@copula@name, interval = FALSE, u = u) : initOpt: DMLE for Gumbel = 0.729065443212314 < 1; is set to 1 Warning message: In bbmle::mle2(minuslogl = nLL, optimizer = "optimize", lower = interval[1], : couldn't invert Hessian > > > ## A simple example with *negative* correlation; > ## Want *more helpful* error message here: > set.seed(7) > u1 <- seq(0,1, by=1/128)[2:127]; u2 <- -u1 + round(rnorm(u1)/4,2); u <- pobs(cbind(u1,u2)) > plot(u) > msg <- tryCatch(fitCopula(gumbelCopula(), data = u), error=function(e)e$message) Warning message: In iTauGumbelCopula(copula, tau) : For the Gumbel copula, tau must be >= 0. Replacing negative values by 0. > ## check error message __FIXME__ want "negative correlation not possible" > ## or NO ERROR and a best fit to tau=0 [and the same for other Archimedean families!] > msg Call: fitCopula(gumbelCopula(), data = u) Fit based on "maximum pseudo-likelihood" and 126 2-dimensional observations. Copula: gumbelCopula alpha 1 The maximized loglikelihood is -1.783e-06 Optimization converged > showProc.time() Time (user system elapsed): 0.64 0 0.64 > > > ## Date: Sat, 14 Nov 2015 20:21:27 -0500 > ## Subject: Replace makePosDef() by nearPD()? > ## From: Marius Hofert > ## To: Ivan Kojadinovic , Jun Yan > ## , Martin Maechler > > ## I currently face the problem of estimating a t-copula to a (504, 413) > ## data matrix. A MWE is: > > require(copula) > set.seed(271) > n <- 504 > d <- 413 > d <- 50 # more realistic -- should become faster !! > U <- matrix(runif(n*d), nrow=n, ncol=d) > if(FALSE) ## takes ??? (> one hour); *control = octr : some trace, not too accurate + system.time(fit <- fitCopula(ellipCopula("t", dim=d, dispstr="un"), + data=U, method="mpl", optim.control = octr)) > > ## > Warning in makePosDef(sigma, delta = 0.001) : > ## > Estimate is not positive-definite. Correction applied. > ## > Process R killed: 9 at Sat Nov 14 19:48:55 2015 > 1 [1] 1 > ## ... so R crashes (this also happens if the data is much 'closer' to a > ## t-copula than my data above, so that's not the problem). I'm wondering > ## about the following. > showProc.time() Time (user system elapsed): 0 0 0 > > ## 4) Implement: first estimating (via pairwise inversion of > ## Kendall's tau) the dispersion matrix and then estimating the d.o.f. > ## via MLE (based on the fitted dispersion matrix). Already in use in > ## vignettes --- would be good to have this in the > ## package as. What we could do is to extend > ## fitCopula.icor(, method="kendall"): if the copula has a d.o.f. parameter, then don't > ## treat it as fixed but estimated it via MLE... (or to get this behavior > ## via a method argument or so?) > > ## 6) After 4), lambda() returns a vector of length 2* d(d-1)/2 ... ok > ## in the bivariate case but not in higher dimensions => want a list > > ###----------- xvCopula() [ <--> ../man/xvCopula.Rd ] --------- > set.seed(12) > x <- rCopula(64, claytonCopula(pi))# "small" n ==> fast > fG <- fitCopula(gumbelCopula(), x) > (v <- c(logL = logLik(fG), + xv.loo = xvCopula(gumbelCopula(), x ), # l.o.o. CV + xv.5 = xvCopula(gumbelCopula(), x, k = 5)))# 5-fold CV logL xv.loo xv.5 32.78368 32.83574 32.24746 > stopifnot(all.equal(unname(v), c(32.783677, 32.835744, 32.247463), + tolerance = 1e-7)) > > ## now also with rotCopula: -- > ## NB: upper=NA / upper=Inf is *better* for optim(*, "L-BFGS-B") than upper= "Large" !!! > ## --- (this is a puzzle in optim() itself !) > if(FALSE)# has been useful: + trace(copula:::fitCopula.ml, + trace = quote({cat("in fitCopula.ml(): copula=\n");print(copula)}) + ) > if(FALSE) ## explicit upper=Inf "works" (*not* finding very good value: + print(ff <- fitCopula(rotCopula(claytonCopula()), upper=NA, x, traceOpt=TRUE)) > if(FALSE) ## explicit upper=Inf or =NA work; upper=100 "fails" giving bad result ! + print(f1c <- fitCopula(rotCopula(claytonCopula()), upper=100, x, traceOpt=TRUE)) > if(FALSE) ## "Brent" instead of default "L-BFGS-B" -- "works" very nicely with upper=100, + print(fB <- fitCopula(rotCopula(claytonCopula()), upper=100, x, + optim.method="Brent",traceOpt=TRUE)) > ## but *not* with upper=Inf or NA (--> error: .. must be finite..) or upper = 1e300 !! > try( fitCopula(rotCopula(claytonCopula()), upper=Inf, x, optim.method="Brent")) Error in optim(start, logL, lower = lower, upper = upper, method = optim.method, : 'lower' and 'upper' must be finite values > try( fitCopula(rotCopula(claytonCopula()), upper=NA, x, optim.method="Brent")) Error in optim(start, logL, lower = lower, upper = upper, method = optim.method, : 'lower' and 'upper' must be finite values > if(FALSE) ## very bad ('start' is not used with "Brent", but the upper bound is!) + print(fBL <- fitCopula(rotCopula(claytonCopula()), upper=1e300, x, + optim.method="Brent",traceOpt=TRUE)) > > ## this "works" with 'Inf', but not with default --> upper = ! > (xvR <- xvCopula(rotCopula(claytonCopula()), upper = Inf, x, k = 8, verbose=FALSE)) [1] 21.05057 > stopifnot(all.equal(xvR, 21.050569, tolerance = 1e-6)) > > ## Now 'copula2 = indepCopula(..) gets 'dim = 3' from the first: > kh.C <- khoudrajiCopula(claytonCopula(2.5, dim=3), shapes = (1:3)/4) > kh.C Khoudraji copula, dim. d = 3, constructed from Clayton copula Independence copula Dimension: 3 Parameters: c1.alpha = 2.50 shape1 = 0.25 shape2 = 0.50 shape3 = 0.75 > getTheta(kh.C, attr=TRUE) # bounds look good c1.alpha shape1 shape2 shape3 2.50 0.25 0.50 0.75 attr(,"param.lowbnd") [1] 0 0 0 0 attr(,"param.upbnd") [1] Inf 1 1 1 > x3 <- cbind(x, x[,1]^2 + runif(nrow(x))) > u3 <- pobs(x3) > if(FALSE) # Error: "non-finite finite-difference [3]" + xvK <- xvCopula(kh.C, u3, k = 8) > if(FALSE) # fails (same Error) + fK <- fitCopula(kh.C, u3) > ## "works" : > p1 <- (1:20)/2 > l1 <- sapply(p1, function(th1) loglikCopula(c(th1, (1:3)/4), u3, kh.C)) > plot(p1, l1, type = "b") # nice maximum at around 6.xx > showProc.time() Time (user system elapsed): 2.16 0.01 2.17 > > ## works two: > fixedParam(kh.C) <- c(FALSE, FALSE, TRUE,TRUE) > getTheta(kh.C) [1] 2.50 0.25 > summary(fK12 <- fitCopula(kh.C, u3)) Call: fitCopula(kh.C, data = u3) Fit based on "maximum pseudo-likelihood" and 64 3-dimensional observations. Khoudraji copula, dim. d = 3, constructed from Clayton copula Independence copula Estimate Std. Error c1.alpha 10.0987 3.855 shape1 0.4647 0.036 The maximized loglikelihood is 36.36 Optimization converged Number of loglikelihood evaluations: function gradient 32 13 Warning message: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" > fixedParam(kh.C) <- FALSE # all free now > fK4 <- fitCopula(kh.C, u3, start = c(coef(fK12), 0.3, 0.5), + optim.method = "L-BFGS-B", traceOpt=TRUE) param= 10.098650142, 0.464681735, 0.300000000, 0.500000000 => logL= 23.089578 param= 10.099650142, 0.464681735, 0.300000000, 0.500000000 => logL= 23.088134 param= 10.097650142, 0.464681735, 0.300000000, 0.500000000 => logL= 23.091021 param= 10.098650142, 0.465681735, 0.300000000, 0.500000000 => logL= 23.01945 param= 10.098650142, 0.463681735, 0.300000000, 0.500000000 => logL= 23.159742 param= 10.098650142, 0.464681735, 0.301000000, 0.500000000 => logL= 23.151984 param= 10.098650142, 0.464681735, 0.299000000, 0.500000000 => logL= 23.026787 param= 10.098650142, 0.464681735, 0.300000000, 0.501000000 => logL= 23.126442 param= 10.098650142, 0.464681735, 0.300000000, 0.499000000 => logL= 23.05218 param= 8.655207058, 0.000000000, 0.999999985, 0.999999985 => logL=1.2873842e-06 param= 8.656207058, 0.000000000, 0.999999985, 0.999999985 => logL=1.2873964e-06 param= 8.654207058, 0.000000000, 0.999999985, 0.999999985 => logL=1.287372e-06 param= 8.655207058, 0.001000000, 0.999999985, 0.999999985 => logL=1.2872785e-06 param= 8.655207058, 0.000000000, 0.999999985, 0.999999985 => logL=1.2873842e-06 param= 8.655207058, 0.000000000, 0.999999985, 0.999999985 => logL=1.2873842e-06 param= 8.655207058, 0.000000000, 0.998999985, 0.999999985 => logL= 0.043308578 param= 8.655207058, 0.000000000, 0.999999985, 0.999999985 => logL=1.2873842e-06 param= 8.655207058, 0.000000000, 0.999999985, 0.998999985 => logL= 0.042972011 param= 9.694071920, 0.334437521, 0.496200843, 0.640143458 => logL= 30.935018 param= 9.695071920, 0.334437521, 0.496200843, 0.640143458 => logL= 30.933812 param= 9.693071920, 0.334437521, 0.496200843, 0.640143458 => logL= 30.936224 param= 9.694071920, 0.335437521, 0.496200843, 0.640143458 => logL= 30.958217 param= 9.694071920, 0.333437521, 0.496200843, 0.640143458 => logL= 30.911813 param= 9.694071920, 0.334437521, 0.497200843, 0.640143458 => logL= 30.860077 param= 9.694071920, 0.334437521, 0.495200843, 0.640143458 => logL= 31.010055 param= 9.694071920, 0.334437521, 0.496200843, 0.641143458 => logL= 30.98665 param= 9.694071920, 0.334437521, 0.496200843, 0.639143458 => logL= 30.88301 param= 9.717282322, 0.339973054, 0.433329322, 0.695425931 => logL= 38.143942 param= 9.718282322, 0.339973054, 0.433329322, 0.695425931 => logL= 38.143093 param= 9.716282322, 0.339973054, 0.433329322, 0.695425931 => logL= 38.14479 param= 9.717282322, 0.340973054, 0.433329322, 0.695425931 => logL= 38.172551 param= 9.717282322, 0.338973054, 0.433329322, 0.695425931 => logL= 38.115019 param= 9.717282322, 0.339973054, 0.434329322, 0.695425931 => logL= 38.065182 param= 9.717282322, 0.339973054, 0.432329322, 0.695425931 => logL= 38.222659 param= 9.717282322, 0.339973054, 0.433329322, 0.696425931 => logL= 38.172548 param= 9.717282322, 0.339973054, 0.433329322, 0.694425931 => logL= 38.114654 param= 9.763552569, 0.486916924, 0.000000000, 0.856683237 => logL= 26.022767 param= 9.764552569, 0.486916924, 0.000000000, 0.856683237 => logL= 26.02278 param= 9.762552569, 0.486916924, 0.000000000, 0.856683237 => logL= 26.022753 param= 9.763552569, 0.487916924, 0.000000000, 0.856683237 => logL= 25.976084 param= 9.763552569, 0.485916924, 0.000000000, 0.856683237 => logL= 26.069472 param= 9.763552569, 0.486916924, 0.001000000, 0.856683237 => logL= 26.039053 param= 9.763552569, 0.486916924, 0.000000000, 0.856683237 => logL= 26.022767 param= 9.763552569, 0.486916924, 0.000000000, 0.857683237 => logL= 26.004989 param= 9.763552569, 0.486916924, 0.000000000, 0.855683237 => logL= 26.040122 param= 9.728732755, 0.376337047, 0.326093924, 0.735332050 => logL= 43.19848 param= 9.729732755, 0.376337047, 0.326093924, 0.735332050 => logL= 43.197943 param= 9.727732755, 0.376337047, 0.326093924, 0.735332050 => logL= 43.199016 param= 9.728732755, 0.377337047, 0.326093924, 0.735332050 => logL= 43.121081 param= 9.728732755, 0.375337047, 0.326093924, 0.735332050 => logL= 43.275322 param= 9.728732755, 0.376337047, 0.327093924, 0.735332050 => logL= 43.222571 param= 9.728732755, 0.376337047, 0.325093924, 0.735332050 => logL= 43.173934 param= 9.728732755, 0.376337047, 0.326093924, 0.736332050 => logL= 43.22142 param= 9.728732755, 0.376337047, 0.326093924, 0.734332050 => logL= 43.174723 param= 9.749442986, 0.325952200, 0.320114380, 0.777971296 => logL= 46.058052 param= 9.750442986, 0.325952200, 0.320114380, 0.777971296 => logL= 46.057544 param= 9.748442986, 0.325952200, 0.320114380, 0.777971296 => logL= 46.058561 param= 9.749442986, 0.326952200, 0.320114380, 0.777971296 => logL= 46.035913 param= 9.749442986, 0.324952200, 0.320114380, 0.777971296 => logL= 46.078644 param= 9.749442986, 0.325952200, 0.321114380, 0.777971296 => logL= 46.037681 param= 9.749442986, 0.325952200, 0.319114380, 0.777971296 => logL= 46.076978 param= 9.749442986, 0.325952200, 0.320114380, 0.778971296 => logL= 46.051501 param= 9.749442986, 0.325952200, 0.320114380, 0.776971296 => logL= 46.064221 param= 9.785681849, 0.291861576, 0.283688457, 0.798643009 => logL= 47.113248 param= 9.786681849, 0.291861576, 0.283688457, 0.798643009 => logL= 47.112515 param= 9.784681849, 0.291861576, 0.283688457, 0.798643009 => logL= 47.11398 param= 9.785681849, 0.292861576, 0.283688457, 0.798643009 => logL= 47.094958 param= 9.785681849, 0.290861576, 0.283688457, 0.798643009 => logL= 47.129986 param= 9.785681849, 0.291861576, 0.284688457, 0.798643009 => logL= 47.097472 param= 9.785681849, 0.291861576, 0.282688457, 0.798643009 => logL= 47.127574 param= 9.785681849, 0.291861576, 0.283688457, 0.799643009 => logL= 47.097792 param= 9.785681849, 0.291861576, 0.283688457, 0.797643009 => logL= 47.128308 param= 9.888749413, 0.186573631, 0.178524589, 0.786413753 => logL= 49.058714 param= 9.889749413, 0.186573631, 0.178524589, 0.786413753 => logL= 49.056789 param= 9.887749413, 0.186573631, 0.178524589, 0.786413753 => logL= 49.060639 param= 9.888749413, 0.187573631, 0.178524589, 0.786413753 => logL= 49.062148 param= 9.888749413, 0.185573631, 0.178524589, 0.786413753 => logL= 49.05363 param= 9.888749413, 0.186573631, 0.179524589, 0.786413753 => logL= 49.052919 param= 9.888749413, 0.186573631, 0.177524589, 0.786413753 => logL= 49.062998 param= 9.888749413, 0.186573631, 0.178524589, 0.787413753 => logL= 49.043471 param= 9.888749413, 0.186573631, 0.178524589, 0.785413753 => logL= 49.073511 param= 9.855863669, 0.184969679, 0.172224207, 0.753028155 => logL= 49.371736 param= 9.856863669, 0.184969679, 0.172224207, 0.753028155 => logL= 49.36966 param= 9.854863669, 0.184969679, 0.172224207, 0.753028155 => logL= 49.373811 param= 9.855863669, 0.185969679, 0.172224207, 0.753028155 => logL= 49.366044 param= 9.855863669, 0.183969679, 0.172224207, 0.753028155 => logL= 49.375839 param= 9.855863669, 0.184969679, 0.173224207, 0.753028155 => logL= 49.374152 param= 9.855863669, 0.184969679, 0.171224207, 0.753028155 => logL= 49.367873 param= 9.855863669, 0.184969679, 0.172224207, 0.754028155 => logL= 49.372295 param= 9.855863669, 0.184969679, 0.172224207, 0.752028155 => logL= 49.370694 param= 9.845107968, 0.184519365, 0.175318702, 0.749193524 => logL= 49.391138 param= 9.846107968, 0.184519365, 0.175318702, 0.749193524 => logL= 49.389056 param= 9.844107968, 0.184519365, 0.175318702, 0.749193524 => logL= 49.393219 param= 9.845107968, 0.185519365, 0.175318702, 0.749193524 => logL= 49.389905 param= 9.845107968, 0.183519365, 0.175318702, 0.749193524 => logL= 49.390738 param= 9.845107968, 0.184519365, 0.176318702, 0.749193524 => logL= 49.388506 param= 9.845107968, 0.184519365, 0.174318702, 0.749193524 => logL= 49.392277 param= 9.845107968, 0.184519365, 0.175318702, 0.750193524 => logL= 49.393296 param= 9.845107968, 0.184519365, 0.175318702, 0.748193524 => logL= 49.388504 param= 9.832908022, 0.184219616, 0.175990314, 0.747375949 => logL= 49.409764 param= 9.833908022, 0.184219616, 0.175990314, 0.747375949 => logL= 49.407676 param= 9.831908022, 0.184219616, 0.175990314, 0.747375949 => logL= 49.411851 param= 9.832908022, 0.185219616, 0.175990314, 0.747375949 => logL= 49.409684 param= 9.832908022, 0.183219616, 0.175990314, 0.747375949 => logL= 49.408202 param= 9.832908022, 0.184219616, 0.176990314, 0.747375949 => logL= 49.405745 param= 9.832908022, 0.184219616, 0.174990314, 0.747375949 => logL= 49.412281 param= 9.832908022, 0.184219616, 0.175990314, 0.748375949 => logL= 49.412662 param= 9.832908022, 0.184219616, 0.175990314, 0.746375949 => logL= 49.406393 param= 9.683271053, 0.180240343, 0.178474265, 0.733512756 => logL= 49.610652 param= 9.684271053, 0.180240343, 0.178474265, 0.733512756 => logL= 49.608496 param= 9.682271053, 0.180240343, 0.178474265, 0.733512756 => logL= 49.612807 param= 9.683271053, 0.181240343, 0.178474265, 0.733512756 => logL= 49.61814 param= 9.683271053, 0.179240343, 0.178474265, 0.733512756 => logL= 49.601505 param= 9.683271053, 0.180240343, 0.179474265, 0.733512756 => logL= 49.597455 param= 9.683271053, 0.180240343, 0.177474265, 0.733512756 => logL= 49.622327 param= 9.683271053, 0.180240343, 0.178474265, 0.734512756 => logL= 49.618592 param= 9.683271053, 0.180240343, 0.178474265, 0.732512756 => logL= 49.602267 param= 9.337910058, 0.170209375, 0.176932799, 0.712201162 => logL= 50.037491 param= 9.338910058, 0.170209375, 0.176932799, 0.712201162 => logL= 50.035185 param= 9.336910058, 0.170209375, 0.176932799, 0.712201162 => logL= 50.039796 param= 9.337910058, 0.171209375, 0.176932799, 0.712201162 => logL= 50.053979 param= 9.337910058, 0.169209375, 0.176932799, 0.712201162 => logL= 50.019415 param= 9.337910058, 0.170209375, 0.177932799, 0.712201162 => logL= 50.013117 param= 9.337910058, 0.170209375, 0.175932799, 0.712201162 => logL= 50.06041 param= 9.337910058, 0.170209375, 0.176932799, 0.713201162 => logL= 50.051693 param= 9.337910058, 0.170209375, 0.176932799, 0.711201162 => logL= 50.022863 param= 8.019821594, 0.131278219, 0.159760853, 0.649339226 => logL= 51.618763 param= 8.020821594, 0.131278219, 0.159760853, 0.649339226 => logL= 51.615995 param= 8.018821594, 0.131278219, 0.159760853, 0.649339226 => logL= 51.62153 param= 8.019821594, 0.132278219, 0.159760853, 0.649339226 => logL= 51.645052 param= 8.019821594, 0.130278219, 0.159760853, 0.649339226 => logL= 51.59122 param= 8.019821594, 0.131278219, 0.160760853, 0.649339226 => logL= 51.573771 param= 8.019821594, 0.131278219, 0.158760853, 0.649339226 => logL= 51.662674 param= 8.019821594, 0.131278219, 0.159760853, 0.650339226 => logL= 51.651408 param= 8.019821594, 0.131278219, 0.159760853, 0.648339226 => logL= 51.58549 param= 4.0497055465, 0.0000000000, 0.0884398043, 0.4972981627 => logL= 57.129611 param= 4.0507055465, 0.0000000000, 0.0884398043, 0.4972981627 => logL= 57.128655 param= 4.0487055465, 0.0000000000, 0.0884398043, 0.4972981627 => logL= 57.130564 param= 4.0497055465, 0.0010000000, 0.0884398043, 0.4972981627 => logL= 57.135645 param= 4.0497055465, 0.0000000000, 0.0884398043, 0.4972981627 => logL= 57.129611 param= 4.0497055465, 0.0000000000, 0.0894398043, 0.4972981627 => logL= 57.061074 param= 4.0497055465, 0.0000000000, 0.0874398043, 0.4972981627 => logL= 57.19778 param= 4.0497055465, 0.0000000000, 0.0884398043, 0.4982981627 => logL= 57.165388 param= 4.0497055465, 0.0000000000, 0.0884398043, 0.4962981627 => logL= 57.093485 param= 4.0180596117, 0.0000000000, 0.0704204478, 0.5051691212 => logL= 58.597989 param= 4.0190596117, 0.0000000000, 0.0704204478, 0.5051691212 => logL= 58.597424 param= 4.0170596117, 0.0000000000, 0.0704204478, 0.5051691212 => logL= 58.598549 param= 4.0180596117, 0.0010000000, 0.0704204478, 0.5051691212 => logL= 58.596053 param= 4.0180596117, 0.0000000000, 0.0704204478, 0.5051691212 => logL= 58.597989 param= 4.0180596117, 0.0000000000, 0.0714204478, 0.5051691212 => logL= 58.536392 param= 4.0180596117, 0.0000000000, 0.0694204478, 0.5051691212 => logL= 58.659142 param= 4.0180596117, 0.0000000000, 0.0704204478, 0.5061691212 => logL= 58.631381 param= 4.0180596117, 0.0000000000, 0.0704204478, 0.5041691212 => logL= 58.564244 param= 5.148262719, 0.000000000, 0.000000000, 0.664442868 => logL= 63.300761 param= 5.149262719, 0.000000000, 0.000000000, 0.664442868 => logL= 63.298988 param= 5.147262719, 0.000000000, 0.000000000, 0.664442868 => logL= 63.302532 param= 5.148262719, 0.001000000, 0.000000000, 0.664442868 => logL= 63.283554 param= 5.148262719, 0.000000000, 0.000000000, 0.664442868 => logL= 63.300761 param= 5.148262719, 0.000000000, 0.001000000, 0.664442868 => logL= 63.28026 param= 5.148262719, 0.000000000, 0.000000000, 0.664442868 => logL= 63.300761 param= 5.148262719, 0.000000000, 0.000000000, 0.665442868 => logL= 63.298913 param= 5.148262719, 0.000000000, 0.000000000, 0.663442868 => logL= 63.30235 param= 4.416664014, 0.000000000, 0.000000000, 0.619546875 => logL= 63.926241 param= 4.417664014, 0.000000000, 0.000000000, 0.619546875 => logL= 63.926177 param= 4.415664014, 0.000000000, 0.000000000, 0.619546875 => logL= 63.926303 param= 4.416664014, 0.001000000, 0.000000000, 0.619546875 => logL= 63.896811 param= 4.416664014, 0.000000000, 0.000000000, 0.619546875 => logL= 63.926241 param= 4.416664014, 0.000000000, 0.001000000, 0.619546875 => logL= 63.901206 param= 4.416664014, 0.000000000, 0.000000000, 0.619546875 => logL= 63.926241 param= 4.416664014, 0.000000000, 0.000000000, 0.620546875 => logL= 63.930667 param= 4.416664014, 0.000000000, 0.000000000, 0.618546875 => logL= 63.921525 param= 4.446598080, 0.000000000, 0.000000000, 0.636720326 => logL= 63.965051 param= 4.447598080, 0.000000000, 0.000000000, 0.636720326 => logL= 63.96505 param= 4.445598080, 0.000000000, 0.000000000, 0.636720326 => logL= 63.965049 param= 4.446598080, 0.001000000, 0.000000000, 0.636720326 => logL= 63.937015 param= 4.446598080, 0.000000000, 0.000000000, 0.636720326 => logL= 63.965051 param= 4.446598080, 0.000000000, 0.001000000, 0.636720326 => logL= 63.939593 param= 4.446598080, 0.000000000, 0.000000000, 0.636720326 => logL= 63.965051 param= 4.446598080, 0.000000000, 0.000000000, 0.637720326 => logL= 63.965058 param= 4.446598080, 0.000000000, 0.000000000, 0.635720326 => logL= 63.964786 param= 4.448405672, 0.000000000, 0.000000000, 0.637274304 => logL= 63.96509 param= 4.449405672, 0.000000000, 0.000000000, 0.637274304 => logL= 63.965088 param= 4.447405672, 0.000000000, 0.000000000, 0.637274304 => logL= 63.965089 param= 4.448405672, 0.001000000, 0.000000000, 0.637274304 => logL= 63.937107 param= 4.448405672, 0.000000000, 0.000000000, 0.637274304 => logL= 63.96509 param= 4.448405672, 0.000000000, 0.001000000, 0.637274304 => logL= 63.939626 param= 4.448405672, 0.000000000, 0.000000000, 0.637274304 => logL= 63.96509 param= 4.448405672, 0.000000000, 0.000000000, 0.638274304 => logL= 63.96497 param= 4.448405672, 0.000000000, 0.000000000, 0.636274304 => logL= 63.964954 param= 4.44824526, 0.00000000, 0.00000000, 0.63730128 => logL= 63.96509 param= 4.44924526, 0.00000000, 0.00000000, 0.63730128 => logL= 63.965089 param= 4.44724526, 0.00000000, 0.00000000, 0.63730128 => logL= 63.965089 param= 4.44824526, 0.00100000, 0.00000000, 0.63730128 => logL= 63.937107 param= 4.44824526, 0.00000000, 0.00000000, 0.63730128 => logL= 63.96509 param= 4.44824526, 0.00000000, 0.00100000, 0.63730128 => logL= 63.939624 param= 4.44824526, 0.00000000, 0.00000000, 0.63730128 => logL= 63.96509 param= 4.44824526, 0.00000000, 0.00000000, 0.63830128 => logL= 63.964962 param= 4.44824526, 0.00000000, 0.00000000, 0.63630128 => logL= 63.964962 param= 4.448253534, 0.000000000, 0.000000000, 0.637301284 => logL= 63.96509 param= 4.449253534, 0.000000000, 0.000000000, 0.637301284 => logL= 63.965089 param= 4.447253534, 0.000000000, 0.000000000, 0.637301284 => logL= 63.965089 param= 4.448253534, 0.001000000, 0.000000000, 0.637301284 => logL= 63.937107 param= 4.448253534, 0.000000000, 0.000000000, 0.637301284 => logL= 63.96509 param= 4.448253534, 0.000000000, 0.001000000, 0.637301284 => logL= 63.939624 param= 4.448253534, 0.000000000, 0.000000000, 0.637301284 => logL= 63.96509 param= 4.448253534, 0.000000000, 0.000000000, 0.638301284 => logL= 63.964962 param= 4.448253534, 0.000000000, 0.000000000, 0.636301284 => logL= 63.964962 > summary(fK4) Call: fitCopula(kh.C, data = u3, ... = pairlist(start = c(coef(fK12), 0.3, 0.5), optim.method = "L-BFGS-B", traceOpt = TRUE)) Fit based on "maximum pseudo-likelihood" and 64 3-dimensional observations. Khoudraji copula, dim. d = 3, constructed from Clayton copula Independence copula Estimate Std. Error c1.alpha 4.4483 1.291 shape1 0.0000 0.129 shape2 0.0000 0.132 shape3 0.6373 0.099 The maximized loglikelihood is 63.97 Optimization converged Number of loglikelihood evaluations: function gradient 23 23 > ## -> shape1 ~= shape2 ~= 0 > ## dput(signif(unname(coef(fK4)), 8)) > c4 <- c(4.4482535, 0, 0, 0.63730128) > all.equal(c4, unname(coef(fK4)), tol=0)# 7.342e-9 [1] "Mean relative difference: 7.346832e-09" > stopifnot(all.equal(c4, unname(coef(fK4)), tol=4e-8)) > showProc.time() Time (user system elapsed): 0.39 0 0.39 > > kh.r.C <- khoudrajiCopula(rotCopula(claytonCopula(1.5)), shapes = c(1,3)/4) > u.krC <- rCopula(999, kh.r.C) > fm.krC <- fitCopula(copula = kh.r.C, data = u.krC, traceOpt=TRUE)# works (w/ warning) start := fitCopStart(..): [1] 1.50 0.25 0.75 param= 1.50, 0.25, 0.75 => logL= 52.909245 param= 1.501, 0.250, 0.750 => logL= 52.903099 param= 1.499, 0.250, 0.750 => logL= 52.915366 param= 1.500, 0.251, 0.750 => logL= 52.915637 param= 1.500, 0.249, 0.750 => logL= 52.902624 param= 1.500, 0.250, 0.751 => logL= 52.927993 param= 1.500, 0.250, 0.749 => logL= 52.88947 param= -4.63349335, 6.75645856, 20.01145957 => logL= -Inf param= 0.27330133, 1.55129171, 4.60229191 => logL= -Inf param= 1.254660266, 0.510258342, 1.520458383 => logL= -Inf param= 1.450932053, 0.302051668, 0.904091677 => logL= 37.44961 param= 1.490186411, 0.260410334, 0.780818335 => logL= 53.089351 param= 1.491186411, 0.260410334, 0.780818335 => logL= 53.086891 param= 1.489186411, 0.260410334, 0.780818335 => logL= 53.09179 param= 1.490186411, 0.261410334, 0.780818335 => logL= 53.096744 param= 1.490186411, 0.259410334, 0.780818335 => logL= 53.081753 param= 1.490186411, 0.260410334, 0.781818335 => logL= 53.072474 param= 1.490186411, 0.260410334, 0.779818335 => logL= 53.104929 param= -3.32896984, 8.52877791, 0.49971412 => logL= -Inf param= 0.526355160, 1.914083848, 0.724597492 => logL= -Inf param= 1.297420160, 0.591145037, 0.769574167 => logL= 44.929319 param= 1.451633161, 0.326557274, 0.778569502 => logL= 53.210324 param= 1.452633161, 0.326557274, 0.778569502 => logL= 53.209483 param= 1.450633161, 0.326557274, 0.778569502 => logL= 53.211142 param= 1.451633161, 0.327557274, 0.778569502 => logL= 53.204031 param= 1.451633161, 0.325557274, 0.778569502 => logL= 53.216432 param= 1.451633161, 0.326557274, 0.779569502 => logL= 53.202006 param= 1.451633161, 0.326557274, 0.777569502 => logL= 53.217475 param= -1.582532907, -0.270853604, 0.436341972 => logL= -Inf param= 0.844799947, 0.207075099, 0.710123996 => logL= 51.216118 param= 1.330266518, 0.302660839, 0.764880400 => logL= 53.458311 param= 1.331266518, 0.302660839, 0.764880400 => logL= 53.459091 param= 1.329266518, 0.302660839, 0.764880400 => logL= 53.457502 param= 1.330266518, 0.303660839, 0.764880400 => logL= 53.453136 param= 1.330266518, 0.301660839, 0.764880400 => logL= 53.463315 param= 1.330266518, 0.302660839, 0.765880400 => logL= 53.449912 param= 1.330266518, 0.302660839, 0.763880400 => logL= 53.465723 param= 1.275969503, 0.259302038, 0.747800771 => logL= 53.629428 param= 1.276969503, 0.259302038, 0.747800771 => logL= 53.629526 param= 1.274969503, 0.259302038, 0.747800771 => logL= 53.629298 param= 1.275969503, 0.260302038, 0.747800771 => logL= 53.629121 param= 1.275969503, 0.258302038, 0.747800771 => logL= 53.62956 param= 1.275969503, 0.259302038, 0.748800771 => logL= 53.627357 param= 1.275969503, 0.259302038, 0.746800771 => logL= 53.630633 param= 1.268669008, 0.255346263, 0.745088913 => logL= 53.631905 param= 1.269669008, 0.255346263, 0.745088913 => logL= 53.631896 param= 1.267669008, 0.255346263, 0.745088913 => logL= 53.631882 param= 1.268669008, 0.256346263, 0.745088913 => logL= 53.631914 param= 1.268669008, 0.254346263, 0.745088913 => logL= 53.631721 param= 1.268669008, 0.255346263, 0.746088913 => logL= 53.631137 param= 1.268669008, 0.255346263, 0.744088913 => logL= 53.631825 param= 1.267188387, 0.255348114, 0.744581699 => logL= 53.631996 param= 1.268188387, 0.255348114, 0.744581699 => logL= 53.631987 param= 1.266188387, 0.255348114, 0.744581699 => logL= 53.631973 param= 1.267188387, 0.256348114, 0.744581699 => logL= 53.631932 param= 1.267188387, 0.254348114, 0.744581699 => logL= 53.631885 param= 1.267188387, 0.255348114, 0.745581699 => logL= 53.631516 param= 1.267188387, 0.255348114, 0.743581699 => logL= 53.631632 param= 1.267440799, 0.255479165, 0.744549145 => logL= 53.632 param= 1.268440799, 0.255479165, 0.744549145 => logL= 53.631981 param= 1.266440799, 0.255479165, 0.744549145 => logL= 53.631985 param= 1.267440799, 0.256479165, 0.744549145 => logL= 53.631915 param= 1.267440799, 0.254479165, 0.744549145 => logL= 53.631908 param= 1.267440799, 0.255479165, 0.745549145 => logL= 53.631582 param= 1.267440799, 0.255479165, 0.743549145 => logL= 53.631573 param= 1.265803838, 0.258598347, 0.749210445 => logL= 53.622226 param= 1.267113407, 0.256103001, 0.745481405 => logL= 53.631617 param= 1.267375321, 0.255603932, 0.744735597 => logL= 53.631985 param= 1.267427704, 0.255504118, 0.744586436 => logL= 53.631999 param= 1.267438180, 0.255484156, 0.744556604 => logL= 53.632 Warning message: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" > xvK.R <- xvCopula(kh.r.C, u.krC, k = 10) Warning messages: 1: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 2: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 3: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 4: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 5: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 6: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 7: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 8: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 9: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" 10: In FUN(copula, u, default, named = TRUE, ...) : no method to get initial parameter for copula of class "khoudrajiExplicitCopula" > showProc.time() Time (user system elapsed): 3.98 0.05 4.04 > > > ## From: Ivan, 27 Jul 2016 08:58 > set.seed(123) > u <- pobs(rCopula(300, joeCopula(4))) > fjc <- fitCopula(joeCopula(), data = u, method = "itau") Warning message: In var.icor(copula, x, method = method) : The variance estimate cannot be computed for a copula of class joeCopula > ## Now a warning. Previously gave 'Error in dCor(cop) :' > ## dTau() method for class "joeCopula" not yet implemented > summary(fjc) Call: fitCopula(joeCopula(), data = u, ... = pairlist(method = "itau")) Fit based on "inversion of Kendall's tau" and 300 2-dimensional observations. Joe copula, dim. d = 2 Estimate Std. Error alpha 3.692 NA > all.equal(target = c(alpha = 4), coef(fjc)) [1] "Mean relative difference: 0.07691788" > fjcM <- fitCopula(joeCopula(), data = u, start = coef(fjc)) > confint(fjcM)# 2.5 % 97.5 % 2.5 % 97.5 % alpha 3.077266 4.294797 > ## alpha 3.077266 4.294797 > summary(fjcM) Call: fitCopula(joeCopula(), data = u, ... = pairlist(start = coef(fjc))) Fit based on "maximum pseudo-likelihood" and 300 2-dimensional observations. Joe copula, dim. d = 2 Estimate Std. Error alpha 3.686 0.311 The maximized loglikelihood is 182.8 Optimization converged Number of loglikelihood evaluations: function gradient 4 4 > if(FALSE) + dput(signif(drop(coef(fjcM, SE=TRUE)), 4)) > stopifnot( + all.equal(c(Estimate = 3.686, `Std. Error` = 0.31), + drop(coef(fjcM, SE=TRUE)), tol = 0.0002) # seen 0.000158 + ) > > > if(!doExtras) q(save="no") ## so the following auto prints > proc.time() user system elapsed 11.18 0.26 11.43