Description Usage Arguments Details Value Author(s) References Examples
Density function, distribution function, random generation,
generator and inverse generator function for the Gumbel Copula with
parameters alpha
. The 4 classic estimation methods for copulas.
1 2 3 4 5 6 7 8 9 | dgumbel(u, v=NULL, alpha, dim=2, warning = FALSE)
pgumbel(u, v=NULL, alpha, dim=2)
rgumbel(n, alpha, dim=2, method=1)
phigumbel(t, alpha=1)
invphigumbel(t, alpha=1)
gumbel.MBE(x, y, marg = "exp")
gumbel.EML(x, y, marg = "exp")
gumbel.IFM(x, y, marg = "exp")
gumbel.CML(x, y)
|
u |
vector of quantiles if argument |
v |
vector of quantiles, needed if |
n |
number of observations. If |
alpha |
parameter of the Copula. Must be greater than 1. |
dim |
an integer specifying the dimension of the copula. |
t |
dummy variable of the generator φ or the inverse generator φ^-1. could be a n-dimensional array. |
method |
an integer code for the method used in simulation. 1 is the common
frailty approach, 2 uses the K function (only valid with |
x,y |
vectors of observations, realizations of random variable X and Y. |
marg |
a character string specifying the marginals of vector (X,Y). It must be either "exp"(default value) or "gamma". |
warning |
a logical (default value |
The Gumbel Hougaard Copula with parameter alpha
is defined by
its generator
φ(t) = (-ln(t))^alpha.
The generator and inverse generator are
implemented in phigumbel
and invphigumbel
respectively.
As an Archimedean copula, its
distribution function is
C(u_1, ...,u_n) = φ^{-1}(φ(u_1)+...+φ(u_n)) = exp(-( (-ln(u_1))^alpha+...+(-ln(u_n))^alpha )^1/α).
pgumbel
and dgumbel
computes the distribution function (expression above) and
the density (n times differentiation of expression above with respect to u_i).
As there is no explicit
formulas for the density of a Gumbel copula, dgumbel
is not yet impemented
for argument dim>3
. This two functions works with a dim
-dimensional array with
the last dimension being equalled to dim
or with a matrix with dim
columns (see examples).
Random generation is carried out with 2 algorithms the common frailty algorithm (method=1) and the 'K' algorithm (method=2). The common frailty algorithm (cf. Marshall & Olkin(1988)) can be sum up in three lines
generate y_1, y_2 from exponential distribution of mean 1,
generate θ from a stable distribution with parameter(1/alpha,1,1,0),
u_1 <- phi(y_1/θ) and u_2 <- phi(y_2/θ).
This algorithm works with any dimension. See Chambers et al(1976) for stable random generation. The 'K' algorithm use the fact the distribution function of random variable C(U,V) is K(t) = t-φ(y)/φ'(t). The algorithm is
generate v_1, t from uniform distribution
generate v_2 from the K distribution i.e. v_2<-K^{-1}(t).
u_1<-φ^{-1}(φ(v_1)v_2) and u_2<-φ^{-1}(φ(v_1)(1-v_2)) .
Warning, the 'K' algorithm does NOT work with dim>2
.
We implements the 4 usual method of estimation for copulas, namely the Exact Maximum
Likelihood (gumbel.EML
), the Inference for Margins (gumbel.IFM
), the
Moment-base Estimation (gumbel.MBE
) and the Canonical Maximum
Likelihood (gumbel.CML
). For the moment, only two types of marginals are
available : exponential distribution (marg="exp"
) and gamma distribution
(marg="gamma"
).
dgumbel
gives the density,
pgumbel
gives the distribution function,
rgumbel
generates random deviates,
phigumbel
gives the generator,
invphigumbel
gives the inverse generator.
gumbel.EML
, gumbel.IFM
, gumbel.MBE
and gumbel.CML
returns the vector of estimates.
Invalid arguments will result in return value NaN
.
A.-L. Caillat, C. Dutang, M. V. Larrieu and T. Nguyen
Nelsen, R. (2006), An Introduction to Copula, Second Edition, Springer.
Marshall & Olkin(1988), Families of multivariate distributions, Journal of the American Statistical Association.
Chambers et al (1976), A method for simulating stable random variables, Journal of the American Statistical Association.
Nelsen, R. (2005), Dependence Modeling with Archimedean Copulas, booklet available at www.lclark.edu/~mathsci/brazil2.pdf.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 | #dim=2
u<-seq(0,1, .1)
v<-u
#classic parametrization
#independance case (alpha=1)
dgumbel(u,v,1)
pgumbel(u,v,1)
#another parametrization
dgumbel(cbind(u,v), alpha=1)
pgumbel(cbind(u,v), alpha=1)
#dim=3 - equivalent parametrization
x <- cbind(u,u,u)
y <- array(u, c(1,11,3))
pgumbel(x, alpha=2, dim=3)
pgumbel(y, alpha=2, dim=3)
dgumbel(x, alpha=2, dim=3)
dgumbel(y, alpha=2, dim=3)
#dim=4
x <- cbind(x,u)
pgumbel(x, alpha=3, dim=4)
y <- array(u, c(2,1,11,4))
pgumbel(y, alpha=3, dim=4)
#independence case
rand <- t(rgumbel(200,1))
plot(rand[1,], rand[2,], col="green", main="Gumbel copula")
#positive dependence
rand <- t(rgumbel(200,2))
plot(rand[1,], rand[2,], col="red", main="Gumbel copula")
#comparison of random generation algorithms
nbsimu <- 10000
#Marshall Olkin algorithm
system.time(rgumbel(nbsimu, 2, dim=2, method=1))[3]
#K algortihm
system.time(rgumbel(nbsimu, 2, dim=2, method=2))[3]
#pseudo animation
## Not run:
anim <-function(n, max=50)
{
for(i in seq(1,max,length.out=n))
{
u <- t(rgumbel(10000, i, method=2))
plot(u[1,], u[2,], col="green", main="Gumbel copula",
xlim=c(0,1), ylim=c(0,1), pch=".")
cat()
}
}
anim(20, 20)
## End(Not run)
#3D plots
#plot the density
x <- seq(.05, .95, length = 30)
y <- x
z <- outer(x, y, dgumbel, alpha=2)
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 100, shade = 0.25, ticktype = "detailed",
xlab = "x", ylab = "y", zlab = "density")
#with wonderful colors
#code of P. Soutiras
zlim <- c(0, max(z))
ncol <- 100
nrz <- nrow(z)
ncz <- ncol(z)
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF",
"cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
couleurs <- tail(jet.colors(1.2*ncol),ncol)
fcol <- couleurs[trunc(z/zlim[2]*(ncol-1))+1]
dim(fcol) <- c(nrz,ncz)
fcol <- fcol[-nrz,-ncz]
persp(x, y, z, col=fcol, zlim=zlim, theta=30, phi=30, ticktype = "detailed",
box = TRUE, xlab = "x", ylab = "y", zlab = "density")
#plot the distribution function
z <- outer(x, y, pgumbel, alpha=2)
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 100, shade = 0.25, ticktype = "detailed",
xlab = "u", ylab = "v", zlab = "cdf")
#parameter estimation
#true value : lambdaX=lambdaY=1, alpha=2
simu <- qexp(rgumbel(200, 2))
gumbel.MBE(simu[,1], simu[,2])
gumbel.IFM(simu[,1], simu[,2])
gumbel.EML(simu[,1], simu[,2])
gumbel.CML(simu[,1], simu[,2])
#true value : lambdaX=lambdaY=1, alphaX=alphaY=2, alpha=3
simu <- qgamma(rgumbel(200, 3), 2, 1)
gumbel.MBE(simu[,1], simu[,2], "gamma")
gumbel.IFM(simu[,1], simu[,2], "gamma")
gumbel.EML(simu[,1], simu[,2], "gamma")
gumbel.CML(simu[,1], simu[,2])
|
[1] NaN 1 1 1 1 1 1 1 1 1 NaN
[1] 0.00 0.01 0.04 0.09 0.16 0.25 0.36 0.49 0.64 0.81 1.00
[1] NaN 1 1 1 1 1 1 1 1 1 NaN
[1] 0.00 0.01 0.04 0.09 0.16 0.25 0.36 0.49 0.64 0.81 1.00
[1] 0.00000000 0.01853315 0.06156706 0.12426461 0.20452561 0.30102374
[7] 0.41280666 0.53914047 0.67943346 0.83319317 1.00000000
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0 0.01853315 0.06156706 0.1242646 0.2045256 0.3010237 0.4128067
[,8] [,9] [,10] [,11]
[1,] 0.5391405 0.6794335 0.8331932 1
[1] NaN 6.922376 3.646759 2.770993 2.510089 2.586175 3.024395
[8] 4.149303 7.366620 23.650259 NaN
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] NaN 6.922376 3.646759 2.770993 2.510089 2.586175 3.024395 4.149303
[,9] [,10] [,11]
[1,] 7.36662 23.65026 NaN
[1] 0.00000000 0.02585824 0.07770595 0.14790462 0.23351222 0.33277038
[7] 0.44446448 0.56768637 0.70172176 0.84598860 1.00000000
, , 1
[,1]
[1,] 0.00000000
[2,] 0.02585824
, , 2
[,1]
[1,] 0.07770595
[2,] 0.14790462
, , 3
[,1]
[1,] 0.2335122
[2,] 0.3327704
, , 4
[,1]
[1,] 0.4444645
[2,] 0.5676864
, , 5
[,1]
[1,] 0.7017218
[2,] 0.8459886
, , 6
[,1]
[1,] 1
[2,] 0
, , 7
[,1]
[1,] 0.02585824
[2,] 0.07770595
, , 8
[,1]
[1,] 0.1479046
[2,] 0.2335122
, , 9
[,1]
[1,] 0.3327704
[2,] 0.4444645
, , 10
[,1]
[1,] 0.5676864
[2,] 0.7017218
, , 11
[,1]
[1,] 0.8459886
[2,] 1.0000000
elapsed
0.006
elapsed
0.199
[1] 1.103536 1.090272 1.757019
[1] 1.109072 1.095742 1.712734
[1] 1.089771 1.102075 1.718454
[1] 1.717899
[1] 1.0106529 2.0246004 0.9562616 1.9088693 3.2273759
[1] 0.8930957 1.7890478 0.8949991 1.7866356 3.2029018
[1] 0.8936754 1.8034034 0.8460049 1.7022972 3.2601631
[1] 3.087149
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.