R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## Copyright (C) 2016 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan
> ##
> ## This program is free software; you can redistribute it and/or modify it under
> ## the terms of the GNU General Public License as published by the Free Software
> ## Foundation; either version 3 of the License, or (at your option) any later
> ## version.
> ##
> ## This program is distributed in the hope that it will be useful, but WITHOUT
> ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
> ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
> ## details.
> ##
> ## You should have received a copy of the GNU General Public License along with
> ## this program; if not, see <http://www.gnu.org/licenses/>.
> 
> ## Khoudraji Copula
> 
> library(copula)
> isExplicit <- copula:::isExplicit
> (doExtras <- copula:::doExtras())
[1] FALSE
> 
> ## if(!dev.interactive(orNone=TRUE)) pdf("khoudraji-ex.pdf")
> 
> ################################################################################
> ## An explicit Khoudraji copula constructed from claytonCopula and indepCoula
> ## with one fixed shape parameter
> kcf <- khoudrajiCopula(copula2 = claytonCopula(6),
+                        shapes = fixParam(c(0.4, 0.95), c(FALSE, TRUE)))
> kcf
Khoudraji copula, dim. d = 2, constructed from
 Independence copula
 Clayton copula 
Dimension:  2 
Parameters (partly fixed, see ':='):
  c2.alpha   = 6.00
  shape1     = 0.40
  shape2    := 0.95
> stopifnot(
+     all.equal(getTheta(kcf, freeOnly = FALSE, named = TRUE) -> th.,
+               c(c2.alpha = 6, shape1 = 0.4, shape2 = 0.95)),
+     identical(getTheta(kcf, named = TRUE), th.[1:2]))
> 
> ## kcf@exprdist$cdf
> stopifnot(isExplicit(kcf))
> 
> set.seed(123)
> u <- rCopula(20, kcf)
> 
> ## cdf/pdf from explicit expression versus from algorithmic implementation
> max(abs(pCopula(u, kcf) - copula:::pKhoudrajiCopula(u, kcf)))
[1] 0
> max(abs(dCopula(u, kcf) - copula:::dKhoudrajiBivCopula(u, kcf)))
[1] 3.375078e-14
> 
> ################################################################################
> ## A nested Khoudraji copula from kcf and a gumbelCopula, still explicit
> k_kcf_g <- khoudrajiCopula(kcf, gumbelCopula(2),
+                            shapes = fixParam(c(.2, .9), c(TRUE, TRUE)))
> k_kcf_g
Khoudraji copula, dim. d = 2, constructed from
 Khoudraji copula constructed from
   Independence copula
   Clayton copula
 Gumbel copula 
Dimension:  2 
Parameters (partly fixed, see ':='):
  c1.c2.alpha   = 6.00
  c1.shape1     = 0.40
  c1.shape2    := 0.95
  c2.alpha      = 2.00
  shape1       := 0.20
  shape2       := 0.90
> stopifnot(isExplicit(k_kcf_g))
> getTheta(k_kcf_g, freeOnly = FALSE, named = TRUE)
c1.c2.alpha   c1.shape1   c1.shape2    c2.alpha      shape1      shape2 
       6.00        0.40        0.95        2.00        0.20        0.90 
> getTheta(k_kcf_g, named = TRUE)
c1.c2.alpha   c1.shape1    c2.alpha 
        6.0         0.4         2.0 
> 
> ## k_kcf_g@exprdist$cdf
> set.seed(456)
> U <- rCopula(100000, k_kcf_g)
> 
> ## check cdf/pdf with C.n and kde2d
> require(MASS)
Loading required package: MASS
> u <- as.matrix(expand.grid((1:9)/10, (1:9)/10))
> ## cdf versus C.n
> max(abs(pCopula(u, k_kcf_g) - C.n(u, U))) ## < 0.0001)
[1] 0.0009825071
> ## pdf versus kde2d
> kde <- kde2d(U[,1], U[,2], n = 9, lims = c(0.1, 0.9, 0.1, 0.9))
> max(abs(dCopula(u, k_kcf_g) / c(kde$z) - 1)) ## relative difference < 0.13
[1] 0.1231915
> 
> ################################################################################
> ## one more nesting: kcf and k_kcf_g
> monster <- khoudrajiCopula(kcf, k_kcf_g,
+                            shapes = fixParam(c(.9, .1), c(TRUE, FALSE)))
> monster; validObject(monster)
Khoudraji copula, dim. d = 2, constructed from
 Khoudraji copula constructed from
   Independence copula
   Clayton copula
 Khoudraji copula constructed from
   Khoudraji copula constructed from
     Independence copula
     Clayton copula
   Gumbel copula 
Dimension:  2 
Parameters (partly fixed, see ':='):
  c1.c2.alpha      = 6.00
  c1.shape1        = 0.40
  c1.shape2       := 0.95
  c2.c1.c2.alpha   = 6.00
  c2.c1.shape1     = 0.40
  c2.c1.shape2    := 0.95
  c2.c2.alpha      = 2.00
  c2.shape1       := 0.20
  c2.shape2       := 0.90
  shape1          := 0.90
  shape2           = 0.10
[1] TRUE
> stopifnot(isExplicit(monster))
> 
> (th.mA <- getTheta(monster, freeOnly = FALSE, named = TRUE))
   c1.c2.alpha      c1.shape1      c1.shape2 c2.c1.c2.alpha   c2.c1.shape1 
          6.00           0.40           0.95           6.00           0.40 
  c2.c1.shape2    c2.c2.alpha      c2.shape1      c2.shape2         shape1 
          0.95           2.00           0.20           0.90           0.90 
        shape2 
          0.10 
> (th.m <- getTheta(monster, named = TRUE))
   c1.c2.alpha      c1.shape1 c2.c1.c2.alpha   c2.c1.shape1    c2.c2.alpha 
           6.0            0.4            6.0            0.4            2.0 
        shape2 
           0.1 
> stopifnot(identical(th.m, th.mA[isFree(monster)]))
> 
> set.seed(488)
> U <- rCopula(100000, monster)
> ## cdf versus C.n
> max(abs(pCopula(u, monster) - C.n(u, U))) ## < 0.002)
[1] 0.001314904
> ## pdf versus kde2d
> kde <- kde2d(U[,1], U[,2], n = 9, lims = c(0.1, 0.9, 0.1, 0.9))
> max(abs(dCopula(u, monster) / c(kde$z) - 1)) ## relative difference < 0.09
[1] 0.08216382
> 
> ################################################################################
> ## khoudrajiExplicitCopula with dim 3
> kcd3 <- khoudrajiCopula(copula1 = indepCopula(dim=3),
+                         copula2 = claytonCopula(6, dim=3),
+                         shapes = c(0.4, 0.95, 0.95))
> kcd3
Khoudraji copula, dim. d = 3, constructed from
 Independence copula
 Clayton copula 
Dimension:  3 
Parameters:
  c2.alpha   = 6.00
  shape1     = 0.40
  shape2     = 0.95
  shape3     = 0.95
> stopifnot(isExplicit(kcd3),
+           inherits(kcd3, "khoudrajiExplicitCopula"))
> 
> set.seed(1712)
> U <- rCopula(100000, kcd3)
> u <- matrix(runif(15), 5, 3)
> ## cdf versus C.n
> max(abs(pCopula(u, kcd3) - C.n(u, U))) # < 0.00007
[1] 0.0006881849
> ## don't know how to check pdf
> (f.v <- dCopula(u, kcd3))
[1] 1.1116773 1.6661734 0.1719080 0.4981574 9.8964415
> stopifnot(
+     all.equal(f.v,
+               c(1.1116773, 1.6661734, 0.1719080, 0.4981574, 9.8964415),
+               tol = 1e-7)
+ )
> 
> kcg <- khoudrajiCopula(claytonCopula(2, dim=3), gumbelCopula(3, dim=3), c(0.2,  0.2, .8))
> kcg
Khoudraji copula, dim. d = 3, constructed from
 Clayton copula
 Gumbel copula 
Dimension:  3 
Parameters:
  c1.alpha   = 2.0
  c2.alpha   = 3.0
  shape1     = 0.2
  shape2     = 0.2
  shape3     = 0.8
> u <- rCopula(10, kcg)
> dCopula(u, kcg)
 [1]  2.780102  2.723464  1.041935  2.027283  1.675718  2.307609  1.203909
 [8] 10.459996  0.386535  3.480239
> 
> 
> kkcgg <- khoudrajiCopula(kcg, gumbelCopula(2, dim = 3), c(.3, .3, .9))
> kkcgg
Khoudraji copula, dim. d = 3, constructed from
 Khoudraji copula constructed from
   Clayton copula
   Gumbel copula
 Gumbel copula 
Dimension:  3 
Parameters:
  c1.c1.alpha   = 2.0
  c1.c2.alpha   = 3.0
  c1.shape1     = 0.2
  c1.shape2     = 0.2
  c1.shape3     = 0.8
  c2.alpha      = 2.0
  shape1        = 0.3
  shape2        = 0.3
  shape3        = 0.9
> u <- rCopula(10, kkcgg)
> dCopula(u, kkcgg)
 [1] 0.3665809 1.5981348 1.0638901 1.9684222 1.5384452 1.5214729 1.3522174
 [8] 0.9529139 1.2088327 0.4565743
> 
> 
> kcgcd3 <- khoudrajiCopula(kcd3, kkcgg, c(.4, .4, .5))
> kcgcd3
Khoudraji copula, dim. d = 3, constructed from
 Khoudraji copula constructed from
   Independence copula
   Clayton copula
 Khoudraji copula constructed from
   Khoudraji copula constructed from
     Clayton copula
     Gumbel copula
   Gumbel copula 
Dimension:  3 
Parameters:
  c1.c2.alpha      = 6.00
  c1.shape1        = 0.40
  c1.shape2        = 0.95
  c1.shape3        = 0.95
  c2.c1.c1.alpha   = 2.00
  c2.c1.c2.alpha   = 3.00
  c2.c1.shape1     = 0.20
  c2.c1.shape2     = 0.20
  c2.c1.shape3     = 0.80
  c2.c2.alpha      = 2.00
  c2.shape1        = 0.30
  c2.shape2        = 0.30
  c2.shape3        = 0.90
  shape1           = 0.40
  shape2           = 0.40
  shape3           = 0.50
> u <- rCopula(10, kcgcd3)
> dCopula(u, kcgcd3)
 [1] 0.5406296 4.7736360 0.5877674 0.7449026 2.2931403 0.3452497 1.5496446
 [8] 0.3569642 2.9840308 1.2362195
> 
> ##--- or use comparederiv() from  ../inst/Rsource/utils.R   via source(..)
> ## dCdu
> all.equal(copula:::dCdu(kcgcd3, u), copula:::dCduNumer(kcgcd3, u, may.warn=FALSE))
[1] "Mean relative difference: 0.0008815003"
> ## dCdtheta
> all.equal(copula:::dCdtheta(kcgcd3, u), copula:::dCdthetaNumer(kcgcd3, u, may.warn=FALSE))
[1] "Mean relative difference: 0.002194312"
> ## dlogcdu
> all.equal(copula:::dlogcdu(kcgcd3, u), copula:::dlogcduNumer(kcgcd3, u, may.warn=FALSE))
[1] TRUE
> ## dlogcdtheta
> all.equal(copula:::dlogcdtheta(kcgcd3, u), copula:::dlogcdthetaNumer(kcgcd3, u, may.warn=FALSE))
[1] TRUE
> 
> #############################################################################
> ## fitting checking
> 
> do1 <- function(n, cop) {
+     U <- pobs(rCopula(n, cop))
+     fit <- try(fitCopula(cop, U, method = "mpl"))
+     p <- nParam(cop, freeOnly=TRUE)
+     if (inherits(fit, "try-error")) rep(NA_real_, 2 * p)
+     else c(coef(fit), sqrt(diag(vcov(fit))))
+ }
> 
> sumsim <- function(sim) {
+     p <- nrow(sim) / 2
+     mm <- apply(sim, 1, mean, na.rm = TRUE)
+     ss <- apply(sim, 1, sd, na.rm = TRUE)
+     data.frame(estimate = mm[1:p], ese = ss[1:p], ase = mm[p + 1:p])
+ }
> 
> if (doExtras) {
+     set.seed(123)
+     nrep <- 50
+     n <- 1000
+ 
+     sim <- replicate(nrep, do1(n, kcg))
+     sumsim(sim)
+     cbind(true = getTheta(kcg), sumsim(sim))
+     ## the average se is greater than empirical for the shape parameters; need more investigation
+ }
> 
> proc.time()
   user  system elapsed 
  12.81    0.31   13.06