# R/BiCopParamDistLp.R In MMDCopula: Robust Estimation of Copulas by Maximum Mean Discrepancy

#### Documented in BiCopParamDistLp

#' Compute the distance between 2 parametric copulas
#'
#' This function uses the numerical integration procedure
#' \code{cubature::\link[cubature]{hcubature}()} to numerical integrate the distance between
#' the distribution or between the densities of two bivariate copulas.
#'
#' @param family family of the first copula.
#'
#' @param par first parameter of the first copula.
#'
#' @param par_p first parameter of the second copula.
#'
#' @param par2 second parameter of the first copula
#' (only useful for two-parameter families of copulas).
#'
#' @param par2_p second parameter of the first copula
#' (only useful for two-parameter families of copulas).
#'
#' @param family_p family of the second copula.
#'
#' @param p determines the \eqn{L_p} distance that is used.
#'
#' @param type type of the functions considered.
#' Can be \code{"cdf"} for the distance between the two cumulative distribution functions
#' or \code{"pdf"} for the distance between the two probability density functions.
#'
#' @param maxEval maximum number of evaluation of the function to be integrated.
#' If 0, then no maximum limit is given. (Only used if \code{p < Inf}).
#'
#' @param truncVal the distance is computed using the supremum or the integral
#' of the function on \eqn{[truncVal, 1 - truncVal]^2}.
#'
#' @return If \code{p < Inf}, it returns a list of four items
#'   \itemize{
#'     \item \code{distance} the value of the distance
#'     \item \code{integral} the value of the integral, which is
#'      the \eqn{p}-th power of the distance.
#'     \item \code{error} the estimated relative error of the integral
#'     \item \code{returnCode} the integer return code of the C routine
#'      This should be 0 if there is no error.
#'   }
#' If \code{p = Inf}, it returns a list of two items
#'   \itemize{
#'     \item \code{distance} the maximum difference between the two copulas
#'     (respectively, between the two copula densities).
#'     \item \code{u_max} the point at which this difference is attained.
#'  }
#'
#' @examples
#' # Distance between the densities of a Gaussian copula with correlation 0.5
#' # and a Gaussian copula with correlation 0.2
#' BiCopParamDistLp(family = 1, par = 0.5, par_p = 0.2, p = 2, type = "cdf", maxEval = 10)
#' BiCopParamDistLp(family = 1, par = 0.5, par_p = 0.2, p = Inf, type = "cdf")
#'
#' # Distance between the cdf of a Student copula
#' # with correlation 0.5 and 4 degrees of freedom
#' # and a Student copula with the same correlation but 20 degrees of freedom
#' BiCopParamDistLp(family = 2, par = 0.5, par_p = 0.5,
#' par2 = 5, par2_p = 20, p = 2, type = "pdf", maxEval = 10)
#'
#' # Distance between the densities of a Gaussian copula with correlation 0.5
#' # and of a Student copula with correlation 0.5 and 15 degrees of freedom
#' BiCopParamDistLp(family = 1, par = 0.5, par_p = 0.5, par2_p = 15,
#' family_p = 2, p = 2, type = "pdf", maxEval = 10)
#'
#' @export
#'
BiCopParamDistLp <- function(family, par, par_p, par2 = par, par2_p = par_p, family_p = family,
p, type, maxEval = 0, truncVal = 0)
{
if (p < Inf) {
switch(
type,
"cdf" = {
toIntegrate <- function(u) {
return ( matrix(( abs(
VineCopula::BiCopCDF(u1 = u[1,], u2 = u[2,],
family = family, par = par, par2 = par2) -
VineCopula::BiCopCDF(u1 = u[1,], u2 = u[2,],
family = family_p, par = par_p, par2 = par2_p)
))^p, nrow = 1, ncol = length(u[1,])) )
}
},

"pdf" = {
toIntegrate <- function(u) {
return ( matrix(( abs(
VineCopula::BiCopPDF(u1 = u[1,], u2 = u[2,],
family = family, par = par, par2 = par2) -
VineCopula::BiCopPDF(u1 = u[1,], u2 = u[2,],
family = family_p, par = par_p, par2 = par2_p)
))^p, nrow = 1, ncol = length(u[1,])) )
}
})

result = cubature::hcubature(f = toIntegrate,
lower = c(truncVal, truncVal),
upper = c(1 - truncVal, 1 - truncVal),
vectorInterface = TRUE, maxEval = maxEval)
result[["distance"]] = result$integral^(1/p) return (result[c("distance", "integral", "error", "returnCode")]) } else if (p == Inf){ # if (maxEval == 0){maxEval <- 10000} # nGrid = floor(sqrt(maxEval)) # univGrid = seq(0 , nGrid/(nGrid+1), length = nGrid) + 1/(2 * nGrid) # u = expand.grid(univGrid , univGrid) # switch( # type, # "cdf" = { # diff_ = # VineCopula::BiCopCDF(u1 = u[,1], u2 = u[,2], # family = family, par = par, par2 = par2) - # VineCopula::BiCopCDF(u1 = u[,1], u2 = u[,2], # family = family_p, par = par_p, par2 = par2_p) # }, # # "pdf" = { # diff_ = # VineCopula::BiCopPDF(u1 = u[,1], u2 = u[,2], # family = family, par = par, par2 = par2) - # VineCopula::BiCopPDF(u1 = u[,1], u2 = u[,2], # family = family_p, par = par_p, par2 = par2_p) # }) # # result = max(abs(diff_)) # which_ = which.max(abs(diff_)) # u_max = as.numeric(u[which_, ]) switch( type, "cdf" = { fn <- function(u) { return ( - (VineCopula::BiCopCDF(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2) - VineCopula::BiCopCDF(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p) )^2 ) } gr <- function(u) { return ( - 2 * (unlist(VineCopula::BiCopHfunc(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2)) - unlist(VineCopula::BiCopHfunc(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p)) ) * (VineCopula::BiCopCDF(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2) - VineCopula::BiCopCDF(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p) )) } }, "pdf" = { fn <- function(u) { return ( - ( VineCopula::BiCopPDF(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2) - VineCopula::BiCopPDF(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p) )^2) } gr <- function(u) { return ( - 2 * c(VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2, deriv = "u1") - VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p, deriv = "u1") , VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2, deriv = "u2") - VineCopula::BiCopDeriv(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p, deriv = "u2") ) * (VineCopula::BiCopPDF(u1 = u[1], u2 = u[2], family = family, par = par, par2 = par2) - VineCopula::BiCopPDF(u1 = u[1], u2 = u[2], family = family_p, par = par_p, par2 = par2_p) )) } }) result = stats::optim(par = c(0.2, 0.2), fn = fn, gr = gr, lower = truncVal, upper = 1 - truncVal, method = "L-BFGS-B") return (list(distance = -result$value, u_max = result\$par ))
}
}


## Try the MMDCopula package in your browser

Any scripts or data that you put into this service are public.

MMDCopula documentation built on April 25, 2022, 5:06 p.m.