| mixCopula | R Documentation |
A mixture of m copulas of dimension d with weights
w_j, j=1,2,\ldots,m is itself a d-dimensional copula,
with cumulative distribution function
C(x) = \sum_{j=1}^m w_j C_j(x),
and (probability) density function
c(x) = \sum_{j=1}^m w_j c_j(x),
where C_j are the CDFs and c_j are the
densities of the m component copulas,
j=1,2,\ldots,m.
mixCopula(coplist, w = NULL)
coplist |
a |
w |
numeric vector of length |
It easy to see that the tail dependencies lambda() and
Spearman's rank correlation rho() can be computed as
mixture of the individual measures.
an object of class mixCopula
khoudrajiCopula, rotCopula also create new
copula models from existing ones.
mC <- mixCopula(list(gumbelCopula(2.5, dim=3),
claytonCopula(pi, dim=3),
tCopula(0.7, dim=3)),
c(2,2,4)/8)
mC
stopifnot(dim(mC) == 3)
set.seed(17)
uM <- rCopula(600, mC)
splom2(uM, main = "mixCopula( (gumbel, clayton, t-Cop) )")
d.uM <- dCopula(uM, mC)
p.uM <- pCopula(uM, mC)
## mix a Gumbel with a rotated Gumbel (with equal weights 1/2):
mGG <- mixCopula(list(gumbelCopula(2), rotCopula(gumbelCopula(1.5))))
rho(mGG) # 0.57886
lambda(mGG)# both lower and upper tail dependency
loglikCopula(c(2.5, pi, rho.1=0.7, df = 4, w = c(2,2,4)/8),
u = uM, copula = mC)
## define (profiled) log-likelihood function of two arguments (df, rho) :
ll.df <- Vectorize(function(df, rho)
loglikCopula(c(2.5, pi, rho.1=rho, df=df, w = c(2,2,4)/8),
uM, mC))
(df. <- 1/rev(seq(1/8, 1/2, length=21)))# in [2, 8] equidistant in 1/. scale
ll. <- ll.df(df., rho = (rh1 <- 0.7))
plot(df., ll., type = "b", main = "loglikCopula((.,.,rho = 0.7, df, ..), u, <mixCopula>)")
if(!exists("Xtras")) Xtras <- copula:::doExtras() ; cat("Xtras: ", Xtras,"\n")
if (Xtras) withAutoprint({
Rhos <- seq(0.55, 0.70, by = 0.01)
ll.m <- matrix(NA, nrow=length(df.), ncol=length(Rhos))
for(k in seq_along(Rhos)) ll.m[,k] <- ll.df(df., rho = Rhos[k])
tit <- "loglikelihood(<tCop>, true param. for rest)"
persp (df., Rhos, ll.m, phi=30, theta = 50, ticktype="detailed", main = tit)
filled.contour(df., Rhos, ll.m, xlab="df", ylab = "rho", main = tit)
})
## fitCopula() example -----------------------------------------------------
## 1) with "fixed" weights :
(mNt <- mixCopula(list(normalCopula(0.95), tCopula(-0.7)), w = c(1, 2) / 3))
set.seed(1452) ; U <- pobs(rCopula(1000, mNt))
(m1 <- mixCopula(list(normalCopula(), tCopula()), w = mNt@w))
getTheta(m1, freeOnly = TRUE, attr = TRUE)
getTheta(m1, named=TRUE)
isFree(m1) # all of them; --> now fix the weights :
fixedParam(m1) <- fx <- c(FALSE, FALSE, FALSE, TRUE, TRUE)
stopifnot(identical(isFree(m1), !fx))
if(Xtras) withAutoprint({ ## time
system.time( # ~ 16 sec (nb-mm4) :
fit <- fitCopula(m1, start = c(0, 0, 10), data = U)
)
fit
summary(fit) #-> incl 'Std.Error' (which seems small for rho1 !)
})
## 2) with "free" weights (possible since copula 1.0-0):
(mNt2 <- mixCopula(list(normalCopula(0.9), tCopula(-0.8)), w = c(1, 3) / 4))
set.seed(1959) ; U2 <- pobs(rCopula(2000, mNt2))
if(Xtras) withAutoprint({ ## time
m2 <- mixCopula(list(normalCopula(), tCopula()), w = mNt@w)
system.time( # ~ 13.5 sec (lynne) :
f2 <- fitCopula(m2, start = c(0, 0, 10, c(1/2, 1/2)), data = U2)
)
f2
summary(f2) # NA for 'Std. Error' as did *not* estimate.variance
summary(f2, orig=FALSE) # default 'orig=TRUE': w-scale; whereas
coef(f2, orig=FALSE) # 'orig=FALSE' => shows 'l-scale' instead
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.