Nothing
######################################################
## ##
## a symmetric copula with cubic quadratic sections ##
## ##
######################################################
# (see Example 3.16 in: Nelsen, Roger B. (2006): An Introduction to Copulas, second edition, Springer)
validCqsCopula <- function(object) {
if (object@dimension != 2)
return("Only copulas with cubic quadratic sections of dimension 2 are supported.")
param <- object@parameters
upper <- object@param.upbnd
lower <- object@param.lowbnd
if (length(param) != length(upper))
return("Parameter and upper bound have non-equal length")
if (length(param) != length(lower))
return("Parameter and lower bound have non-equal length")
if (any(is.na(param) | param > upper | param < lower))
return("Parameter value out of bound")
if (object@fixed != ""){
if(!("a" %in% object@fixed | "b" %in% object@fixed))
return("The slot fixed may only refer to \"a\" or \"b\".")
if ("a" %in% object@fixed & "b" %in% object@fixed)
return("Only one of the parameters may be kept fixed.")
}
else return (TRUE)
}
setClass("cqsCopula",
representation = representation("copula",fixed="character"),
validity = validCqsCopula,
contains = list("copula")
)
# constructor
cqsCopula <- function (param=c(0,0), fixed="") {
new("cqsCopula", dimension = as.integer(2), parameters = param,
param.names = c("a", "b"), param.lowbnd = c(limA(param[2]),-1),
param.upbnd = c(1, 1),
fullname = "copula family with cubic quadratic sections", fixed=fixed)
}
## printing
setMethod("describeCop", c("cqsCopula", "character"),
function(x, kind = c("short", "very short", "long"), prefix = "", ...) {
kind <- match.arg(kind)
if(kind == "very short") # e.g. for show() which has more parts
return(paste0(prefix, "CQS copula"))
name <- "cubic-quadratic sections"
d <- dim(x)
ch <- paste0(prefix, name, " copula, dim. d = ", d)
switch(kind <- match.arg(kind),
short = ch,
long = paste0(ch, "\n", prefix, " param.: ",
capture.output(str(x@parameters,
give.head=FALSE))),
stop("invalid 'kind': ", kind))
})
## density ##
dCQSec <- function (u, copula, log=F) {
a <- copula@parameters[1]
b <- copula@parameters[2]
if (!is.matrix(u)) u <- matrix(u, ncol = 2)
u1 <- u[, 1]
u2 <- u[, 2]
if (log)
return(log(pmax(1-b*(1-2*u2)*(1-2*u1)+(b-a)*(1-u2)*(1-3*u2)*(1-u1)*(1-3*u1),0)))
else
return(pmax(1-b*(1-2*u2)*(1-2*u1)+(b-a)*(1-u2)*(1-3*u2)*(1-u1)*(1-3*u1),0))
}
setMethod("dCopula", signature("numeric", "cqsCopula"),
function(u, copula, log) {
dCQSec(matrix(u,ncol=copula@dimension), copula, log)
})
setMethod("dCopula", signature("matrix", "cqsCopula"), dCQSec)
## jcdf ##
pCQSec <- function (u, copula) {
a <- copula@parameters[1]
b <- copula@parameters[2]
if (!is.matrix(u))
u <- matrix(u, ncol = 2)
u1 <- u[, 1]
u2 <- u[, 2]
return(u1*u2*(1- b*(1-u1)*(1-u2) + (b-a)*(1-u2)^2*(1-u1)^2))
}
setMethod("pCopula", signature("numeric", "cqsCopula"),
function(u, copula, ...) {
pCQSec(matrix(u,ncol=copula@dimension), copula)
})
setMethod("pCopula", signature("matrix","cqsCopula"), pCQSec)
## partial derivatives ##
# solves a*x^3 + b*x^2 + c*x + d = 0
solveCubicEq <- function(a,b,c,d){
eps <- .Machine$double.eps
# using the reduced equation z^3 + 3 * p * z + q = 0 with:
p <- 3*a*c-b^2
q <- 2*b^3-9*a*b*c+27*a^2*d
D <- q^2+4*p^3
z <- matrix(NA,nrow=length(D),ncol=3)
ind <- abs(D) <= eps
if(any(ind)){
z[ind,1] <- 0.5*(-4*q[ind])^(1/3)
z[ind,2] <- -z[ind,1]
}
ind <- D > eps
if(any(ind)){
cubeRad <- -4*q[ind]+4*sqrt(D[ind])
r1 <- sign(cubeRad)*abs(cubeRad)^(1/3)
cubeRad <- -4*q[ind]-4*sqrt(D[ind])
r2 <- sign(cubeRad)*abs(cubeRad)^(1/3)
z[ind,1] <- 0.5*(r1+r2)
}
ind <- D < eps
if(any(ind)){
phi <- acos(-q[ind]/(2*sqrt(-p[ind]^3)))
triple <- NULL
triple <- sqrt(-p[ind])*2*cos(phi/3)
triple <- cbind(triple,sqrt(-p[ind])*2*cos((phi+2*pi)/3))
triple <- cbind(triple,sqrt(-p[ind])*2*cos((phi+4*pi)/3))
z[ind,] <- triple
}
return((z-b)/(3*a))
}
## partial derivative ddu pCQSec
dduCQSec <- function (u, copula) {
a <- copula@parameters[1]
b <- copula@parameters[2]
u1 <- u[, 1]
u2 <- u[, 2]
return(u2-b*(u2-u2^2-2*u1*u2+2*u1*u2^2)+(b-a)*(u2-4*u1*u2+3*u1^2*u2-2*u2^2+8*u1*u2^2-6*u1^2*u2^2+u2^3-4*u1*u2^3+3*u1^2*u2^3))
}
setMethod("dduCopula", signature("numeric","cqsCopula"),
function(u, copula, ...) {
dduCQSec(matrix(u,ncol=copula@dimension), copula)
})
setMethod("dduCopula", signature("matrix","cqsCopula"), dduCQSec)
## inverse partial derivative ddu
## inverse partial derivative ddu
# seems to be accurate (1.4e-05 is the max out of 1000 random CQSec-copulas for 1000 random pairs (u,v) each.)
invdduCQSec <- function (u, copula, y) {
stopifnot(length(u) == length(y))
a <- copula@parameters[1]
b <- copula@parameters[2]
if (a != b) { # the cubic case
# solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
usq <- u^2
c3 <- (b-a)*(1-4*u+3*usq)
c2 <- (b-a)*(-2+8*u-6*u^2)-b*(-1+2*u)
c1 <- (b-a)*(1-4*u+3*u^2)-b*(1-2*u)+1
c0 <- -y
v <- solveCubicEq(c3,c2,c1,c0)
filter <- function(vec){
vec <- vec[!is.na(vec)]
return(vec[vec >= 0 & vec <= 1])
}
return(apply(v,1,filter))
}
if(a==0) # and b==0 obvioulsy as well: the independent case
return(y)
# the qudratic cases remain
v <- y
uR <- u[u != 0.5]
v[u != 0.5] <- (-sqrt((-2*b*uR+b-1)^2-4*y[u != 0.5]*(2*b*uR-b))+2*b*uR-b+1)/(2*b*(2*uR-1))
return(v)
}
setMethod("invdduCopula", signature("numeric","cqsCopula","numeric"), invdduCQSec)
## partial derivative ddv
ddvCQSec <- function (u, copula) {
a <- copula@parameters[1]
b <- copula@parameters[2]
if (!is.matrix(u)) u <- matrix(u, ncol = 2)
u1 <- u[, 1]
u2 <- u[, 2]
return(u1-b*(u1-2*u1*u2-u1^2+2*u1^2*u2)+(b-a)*(u1-2*u1^2+u1^3-4*u1*u2+8*u1^2*u2-4*u1^3*u2+3*u1*u2^2-6*u1^2*u2^2+3*u1^3*u2^2))
}
setMethod("ddvCopula", signature("numeric","cqsCopula"),
function(u, copula, ...) {
ddvCQSec(matrix(u,ncol=copula@dimension), copula)
})
setMethod("ddvCopula", signature("matrix","cqsCopula"), ddvCQSec)
## inverse partial derivative ddv
invddvCQSec <- function (v, copula, y) {
stopifnot(length(v)==length(y))
a <- copula@parameters[1]
b <- copula@parameters[2]
if (a != b) { # the cubic case
# solving the cubic equation: u^3 * c3 + u^2 * c2 + u * c1 + c0 = 0
vsq <- v^2
c3 <- (b-a)*(1-4*v+3*vsq)
c2 <- (b-a)*(-2+8*v-6*vsq)-b*(-1+2*v)
c1 <- (b-a)*(1-4*v+3*vsq)-b*(1-2*v)+1
c0 <- -y
u <- solveCubicEq(c3,c2,c1,c0)
filter <- function(vec){
vec <- vec[!is.na(vec)]
return(vec[vec >= 0 & vec <= 1])
}
return(apply(u,1,filter))
}
if(a==0) # and b==0 obvioulsy as well: the independent case
return(y)
# the qudratic cases remain
u <- y
vR <- v[v != 0.5]
u[v != 0.5] <- (-sqrt((-2*b*vR+b-1)^2-4*y[v != 0.5]*(2*b*vR-b))+2*b*vR-b+1)/(2*b*(2*vR-1))
return(u)
}
setMethod("invddvCopula", signature("numeric","cqsCopula","numeric"), invddvCQSec)
## random number generator
rCQSec <- function (n, copula) {
u <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
res <- cbind(u, invdduCQSec(u, copula, y))
colnames(res) <- c("u","v")
return(res)
}
setMethod("rCopula", signature("numeric","cqsCopula"), rCQSec)
## fitment
fitCopula.cqs <- function (copula, data, method = "ml", start=c(0,0),
lower=c(-3,-1), upper=c(1,1),
optim.method="L-BFGS-B", optim.control=list(),
estimate.variance = FALSE) {
fit <- switch(method,
ml=fitCQSec.ml(copula, data, start, lower, upper, optim.control, optim.method),
itau=fitCQSec.itau(copula, data, estimate.variance),
irho=fitCQSec.irho(copula, data, estimate.variance),
stop("Implemented methods for copulas in the spcopula package are: ml, itau, and irho."))
return(fit)
}
setMethod("fitCopula", signature("cqsCopula"), fitCopula.cqs)
## Fits the copula with cubic and quadratic sections acoording to a measure of association.
## It performs a maximum likelihood evaluation over all possible pairs of
## parameters a and b generating a copula with the given mesaure of
## association.
#
# moa
# measure of association, according to method
# data
# the bivariate data set as a 2-column matrix within the unitsquare
# method
# one of kendall or spearman according to the calculation of moa
fitCQSec.itau <- function(copula, data, estimate.variance=FALSE, tau=NULL) {
if(is.null(tau))
tau <- TauMatrix(data)[1,2]
if(copula@fixed=="a")
esti <- c(copula@parameters[1], iTauCQSec.a(copula@parameters[1], tau))
if(copula@fixed=="b")
esti <- c(iTauCQSec.b(copula@parameters[2], tau),copula@parameters[2])
else
esti <- fitCQSec.moa(tau, data, method="itau")
copula <- cqsCopula(esti, fixed=copula@fixed)
new("fitCopula", estimate = esti, var.est = matrix(NA),
method = "Inversion of Kendall's tau and MLE",
loglik = sum(dCopula(data, copula, log=T)),
fitting.stats=list(convergence = as.integer(NA)), nsample = nrow(data),
copula=copula)
}
fitCQSec.irho <- function(copula, data, estimate.variance=FALSE, rho=NULL){
if(is.null(rho))
rho <- cor(data,method="spearman")[1,2]
if(copula@fixed=="a")
esti <- c(copula@parameters[1], iRhoCQSec.a(copula@parameters[1],rho))
if(copula@fixed=="b")
esti <- c(iRhoCQSec.b(copula@parameters[2],rho),copula@parameters[2])
else
esti <- fitCQSec.moa(rho, data, method="irho")
copula <- cqsCopula(esti, fixed=copula@fixed)
new("fitCopula", estimate = esti, var.est = matrix(NA),
method = "Inversion of Spearman's rho and MLE",
loglik = sum(dCopula(data, copula, log=T)),
fitting.stats=list(convergence = as.integer(NA)), nsample = nrow(data),
copula=copula)
}
fitCQSec.moa <- function(moa, data, method="itau", tol=.Machine$double.eps^.5) {
smpl <- as.matrix(data)
iRho.CQS <- function(p) {
iRhoCQSec.b(p,moa)
}
iTau.CQS <- function(p) {
iTauCQSec.b(p,moa)
}
iFun <- switch(method, itau=iTau.CQS, irho=iRho.CQS)
sec <- function (parameters) {
res <- NULL
for(param in parameters) {
res <- rbind(res, -sum(dCQSec(smpl, cqsCopula(c(iFun(param),param)),log=T)))
}
return(res)
}
b <- optimize(sec,c(-1,1), tol=tol)$minimum
return(c(iFun(b),b))
}
# maximum log-likelihood estimation of a and b using optim
fitCQSec.ml <- function(copula, data, start, lower, upper, optim.control, optim.method) {
if(length(start)!=2) stop("Start values need to have same length as parameters.")
if (copula@fixed=="") {
optFun <- function(param=c(0,0)) {
if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2]))
return(100)
return(-sum(log( dCQSec(data, cqsCopula(param)) )))
}
optimized <- optim(par=start, fn=optFun, method = optim.method,
lower=lower, upper=upper, control = optim.control)
return(new("fitCopula", estimate = optimized$par, var.est = matrix(NA),
method = "Numerical MLE over the full range.",
loglik = -optimized$value, fitting.stats= optimized,
nsample = nrow(data), copula=cqsCopula(optimized$par)))
} else {
optFunFixed <- function(p) {
param <- switch(copula@fixed, a=c(copula@parameters[1],p),
b=c(p,copula@parameters[2]))
if(any(param > 1) | param[2] < -1 | param[1] < limA(param[2]))
return(100)
return(-sum(log( dCQSec(data, cqsCopula(param)) )))
}
optimized <- optimise(optFunFixed, c(-3,1))
optPar <- switch(copula@fixed, a=c(copula@parameters[1],optimized$minimum),
b=c(optimized$minimum,copula@parameters[2]))
return(new("fitCopula", estimate = optimized$minimum, var.est = matrix(NA),
method = "Numerical MLE over the full range.",
loglik = -optimized$objective, fitting.stats=list(),
nsample = nrow(data), copula=cqsCopula(optPar)))
}
}
####
iTauCQSec.a <- function(a, tau=0) {
limits <- limB(a)
min(max(limits[1],0.5*(sqrt(a^2-250*a-1800*tau+5626)+a-75)),limits[2])
}
iTauCQSec.b <- function(b,tau=0) {
min(max(limA(b),(b^2 + 75*b + 450*tau)/(b - 25)),1)
}
setMethod("iTau",signature="cqsCopula",
function(copula, tau) {
switch(copula@fixed,
a=c(copula@parameters[1],iTauCQSec.a(copula@parameters[1],tau)),
b=c(iTauCQSec.b(copula@parameters[2],tau),copula@parameters[2]),
stop("iTau may only be used for cqsCopula with one parameter fixed."))
})
####
tauCQSec <- function(copula){
a <- copula@parameters[1]
b <- copula@parameters[2]
return( (a*b - 25*a - b^2 - 75*b)/450 )
}
setMethod("tau",signature("cqsCopula"),tauCQSec)
####
# find parameter "a" for parameter "b" under a given measure of association "rho"
# it may return a value exceeding the limit of "a" which may result in an invalid copula.
iRhoCQSec.a <- function(a, rho=0) {
limits <- limB(a)
min(max(limits[1],-a/3 - 4*rho),limits[2])
}
iRhoCQSec.b <- function(b, rho=0) {
min(max(limA(b),-3*b - 12*rho),1)
}
setMethod("iRho",signature="cqsCopula",
function(copula, rho) {
switch(copula@fixed,
a=function(copula, rho) c(copula@parameters[1],iRhoCQSec.a(copula@parameters[1],rho)),
b=function(copula, rho) c(iRhoCQSec.b(copula@parameters[2],rho),copula@parameters[2]),
stop("iRho may only be used for cqsCopula with one parameter fixed."))
})
####
rhoCQSec <- function(copula){
a <- copula@parameters[1]
b <- copula@parameters[2]
return( -(a+3*b)/12 )
}
setMethod("rho",signature("cqsCopula"),rhoCQSec)
setMethod("lambda", signature("cqsCopula"),
function(copula, ...) c(lower = 0, upper = 0))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.