[R] Error while fitting gumbel copula

Tue Aug 23 13:22:57 CEST 2016

Dear Martin,

Thank you for the solutions.

I've tried as per your codes and my problem solved :)

Thanks a lot!

Best regards,
Isaudin

On Tue, Aug 23, 2016 at 9:06 AM, Martin Maechler <maechler at stat.math.ethz.ch
> wrote:

> >>>>> Isaudin Ismail <isaudin at gmail.com>
> >>>>>     on Thu, 18 Aug 2016 17:03:50 +0100 writes:
>
>     > Dear Martin, Following my earlier question on "error while
>     > fitting gumbel copula", I have also crated a new gist at
>     > https://gist.github.com/anonymous/0bb8aba7adee550d40b840a47d8b7e25
>     > for easy checking and copying codes.
>
>     > I got no problem fitting other Archimedean copulas except
>     > gumbel copula as per my code I used above.
>
>     > Appreciate your kind help.
>
>     > Many thanks, Isaudin
>
>
>     > On Mon, Aug 15, 2016 at 4:28 PM, Isaudin Ismail
>     > <isaudin at gmail.com> wrote:
>
>     >> Dear Dr. Martin,
>     >>
>     >> I'm glad that you replied to my queries.
>     >>
>     >> As advised, I have prepared the following:
>
> MM: I'm including (cut'n'pasting) my commented and augmented version here:
>
> ------------------------------------------------------------------
>
> ## From: Isaudin Ismail <isaudin at gmail.com>
> ## To: Martin Maechler <maechler at stat.math.ethz.ch>
> ## CC: <r-help at r-project.org>
> ## Subject: Re: [R] Error while fitting gumbel copula
> ## Date: Mon, 15 Aug 2016 16:28:14 +0100
>
> ## Dear Dr. Martin,
>
> ## I'm glad that you replied to my queries.
> ## As advised, I have prepared the following:
>
> library(copula)
>
> ## 5 series of data, A, B, C, D and E
> A <- c(0.849420849, 0.900652985, 0.97144217, 0.817888428, 0.877901578,
>        1.070040669, 0.889742431, 0.87588968, 0.853541938, 0.848664688,
>        0.876830319, 0.749582638, 0.818515498, 0.890997174, 0.794766966,
>        0.784794851, 0.814858959, 1.074396518, 0.83752495, 0.894341116,
>        0.880375293, 0.900816803)
>
> B <- c(0.479850746, 0.652111668, 1.880607815, 0.579902303, 0.50669344,
>        0.747560182, 0.701754386, 0.48969697, 0.346751006, 0.379234973,
>        0.862691466, 0.328280188, 0.317312661, 0.534438115, 0.487002653,
>        0.335043612, 0.373346897, 0.627520161, 0.792114695, 0.938253012,
>        0.444553967, 0.625972763)
>
> C <- c(0.693491124, 0.866523143, 4.585714286, 1.512055109, 0.387755102,
>        0.513435701, 0.76252505, -0.113113113, 0.338521401, 0.333951763,
>        0.668755595, 0.401273885, 0.419868791, 0.272885789, 0.541541542,
>        0.32751938, 0.386409736, 0.957446809, 0.861195542, 1.531632653,
>        0.431610942, 1.226470588)
>
> D <- c(0.807792208, 0.548547718, 0.738232865, 0.542247744, 1.088964927,
>        0.862385321, 0.60720268, 1.000816993, 0.699289661, 0.41723356,
>        0.604037267, 0.605003791, 0.698940998, 0.764792899, 0.647897898,
>        0.825256975, 0.767476085, 0.941391941, 0.889547813, 0.324503311,
>        0.942435424, 0.740686633)
>
> E <- c(1.077598829, 0.318507891, 1.152616279, 0.930397727, 1.515994437,
>        0.940689655, 0.880886427, 1.054274084, 1.067282322, 0.677419355,
>        0.966233766, 0.761029412, 1.05734767, 0.615925059, 1.061988304,
>        1.07184241, 1.058890147, 1.123873874, 1.304891923, -0.069584736,
>        1.172757475, 0.501096491)
>
> require(copula)
> gumbel.copula <- gumbelCopula(dim = 2)
> p <- pobs(cbind(D + E, A + B+ C ))
>
> fit.gumbel <- fitCopula(gumbel.copula, p, method = "ml")
>
> ## The error is here when trying to fit the gumbel copula
>
> # I got the following error:
> ## 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]
>
> ## MM:  my version of copula gives the error message  "some tau < 0"
> ## --                                                      ---------
> ## and indeed:
> (tau.p <- cor(p[,1], p[,2], method="kendall"))
> ## [1] -0.1428571
> ##     ^---------  Kendall's tau is =  - 1/7  < 0  ... and that is not
> good for Gumbel!
>
> plot(p)
>
> ##----------------------------------------------------------
> -------------------
>
> So, you tried fitting to *negatively* correlated data, and if
> you use the default instead of "ml" the copula is fit, and uses
> param = 1 (which corresponds to the *independence* copula:
> Because among all the (weakly) positively correlated gumbel
> copulas, the boundary case, param = 1 (<==> tau = 0) is the best
> fitting.
>
> What you can do is to "rotate" the data (actually mirror it),
> and fit a gumbel copula, which now works nice and easily :
>
> p2 <- p; p2[,2] <- 1-p2[,2]
> (tau.p2 <- cor(p2, method="kendall"))
> ## --> now positively correlated
> ## --->
> gumb.ml.p2 <- fitCopula(gumbel.copula, p2, method = "ml")
> summary(gumb.ml.p2) # looks fine now :
>
>   Call: fitCopula(copula, data = data, method = "ml")
>   Fit based on "maximum likelihood" and 22 2-dimensional observations.
>   Gumbel copula, dim. d = 2
>           Estimate Std. Error
>   param    1.121      0.209
>   The maximized loglikelihood is 0.1839
>   Optimization converged
>   Number of loglikelihood evaluations:
>   function gradient
>          6        6
>
>
> ---
>
> The next version of copula --- or the R-forge one, if you are
> interested will support fitting "rotated" copulas,
> so you fit a rotated gumbel copula, flipping, the 2nd coordinate
> (but not the first):
>
> rotGcop <- rotCopula(gumbelCopula(), flip=c(FALSE,TRUE))
>
> f2 <- fitCopula(rotGcop, data = p) # default method: "mlp"
> summary(f2)
>
>   ## Call: fitCopula(copula, data = data)
>   ## Fit based on "maximum pseudo-likelihood" and 22 2-dimensional
> observations.
>   ## rotCopula copula:  Gumbel copula, dim. d = 2
>   ##       Estimate Std. Error
>   ## param    1.121      0.225
>   ## The maximized loglikelihood is 0.1839
>   ## Optimization converged
>   ## Number of loglikelihood evaluations:
>   ## function gradient
>   ##        6        6
>
> -------
>
> Best regards,
>
> Martin Maechler,
> ETH Zurich
>

	[[alternative HTML version deleted]]