gumbel: The Gumbel Hougaard Copula
In gumbel: The Gumbel-Hougaard Copula
Description
Usage
Arguments
Details
Value
Author(s)
References
Examples
Description
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.
Usage
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)
Arguments
u
vector of quantiles if argument v is provided
or matrix of quantiles if argument v is not provided
v
vector of quantiles, needed if u is not a matrix
n
number of observations. If length(n) > 1, the length is
taken to be the number required.
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 dim=2).
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 FALSE) if you want warnings.
Details
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").
Value
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.
Author(s)
A.-L. Caillat, C. Dutang, M. V. Larrieu and T. Nguyen
References
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.
Examples
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])
Example output
[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
gumbel documentation built on May 1, 2019, 8:04 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Author(s) References Examples
Description
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.
Usage
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)
|
Arguments
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 |
Details
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").
Value
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.
Author(s)
A.-L. Caillat, C. Dutang, M. V. Larrieu and T. Nguyen
References
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.
Examples
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])
|
Example output
[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
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.