R/internal_PIOS_normal.R

In gofCopula: Goodness-of-Fit Tests for Copulae

# Internal functions for the derivation of the normal copula for the PIOSRn and 
# PIOSTn test statistics. Called by .Rn and .Tn, the internal functions for
# their test statistics.

# density
.gadens <- function(u, sig, dims) {
  res.f <- exp(-u %*% (solve(sig) - diag(1, dims)) %*% t(t(u)) / 2) / 
    sqrt(det(sig))
  return(res.f)
}

# density bootstrap
.opt.ga <- function(b, psn.sample, m, dims) {
  sig <- cor(psn.sample[-(((b - 1) * m + 1):(b * m)), ])
  if (m == 1) {
    result <- .gadens(psn.sample[((b - 1) * m + 1):(b * m), ], sig, 
                      dims = dims)
  } else {
    uu <- psn.sample[((b - 1) * m + 1):(b * m), ]
    result <- apply(uu, 1, FUN = .gadens, sig = sig)
  }
  return(result)
}

# first and second derivatives wrt theta of the log-density
.V.ga <- function(u, sig.inv, dims) {
  return(sig.inv %*% (u %*% t(u) %*% sig.inv - diag(1, dims)))
}

.S.ga.12 <- function(u, sig) {
  res.f <- (-1 + sig[1, 2]^4 - 2 * u[1] * u[2] * sig[1, 2] * 
              (3 + sig[1, 2]^2) + (u[1]^2) * (1 + 3 * (sig[1, 2]^2)) + 
              (u[2]^2) * (1 + 3 * (sig[1, 2]^2))) / ((-1 + sig[1, 2]^2)^3)
  return(res.f)
}

.S.ga.12.12 <- function(u, sig) {
  res.f <- (-((-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * 
                 sig[1, 3] * sig[2, 3] + (sig[2, 3]^2)) * (-1 + 3 * 
                  (sig[1, 2]^2) - (u[2]^2) * (-1 + (sig[1, 3]^2)) - 2 * 
                  u[2] * u[3] * sig[2, 3] + (sig[2, 3]^2) + (u[3]^2) * 
                  (sig[2, 3]^2) + 2 * sig[1, 2] * sig[1, 3] * (u[2] * u[3] - 
                  (3 + (u[3]^2)) * sig[2, 3]) - (u[1]^2) * (-1 + 
                  (sig[2, 3]^2)) + (sig[1, 3]^2) * (1 + (u[3]^2) + 2 * 
                  (sig[2, 3]^2)) - 2 * u[1] * (u[3] * (sig[1, 3] - sig[1, 2] * 
                  sig[2, 3]) + u[2] * (sig[1, 2] - sig[1, 3] * sig[2, 3])))) + 
                  4 * (sig[1, 2] - sig[1, 3] * sig[2, 3]) * ((sig[1, 2]^3) + 
                  (sig[1, 2]^2) * sig[1, 3] * (u[2] * u[3] - (3 + (u[3]^2)) * 
                  sig[2, 3]) - (u[1]^2) * (sig[1, 2] - sig[1, 3] * sig[2, 3]) * 
                  (-1 + (sig[2, 3]^2)) + sig[1, 3] * ((u[2]^2) * (-1 + 
                  (sig[1, 3]^2)) * sig[2, 3] + u[2] * u[3] * (1 - 
                  (sig[1, 3]^2) + (sig[2, 3]^2)) - sig[2, 3] * (-1 + (u[3]^2) + 
                  (sig[1, 3]^2) + (sig[2, 3]^2))) + sig[1, 2] * 
                  (-1 - (u[2]^2) * (-1 + (sig[1, 3]^2)) - 2 * u[2] * u[3] * 
                  sig[2, 3] + (sig[2, 3]^2) + (u[3]^2) * (sig[2, 3]^2) + 
                  (sig[1, 3]^2) * (1 + (u[3]^2) + 2 * (sig[2, 3]^2))) + u[1] * 
                  (u[2] * (-1 - (sig[1, 2]^2) + (sig[1, 3]^2) + 2 * sig[1, 2] * 
                  sig[1, 3] * sig[2, 3] + (sig[2, 3]^2) - 2 * (sig[1, 3]^2) * 
                  (sig[2, 3]^2)) + u[3] * (-2 * sig[1, 2] * sig[1, 3] + 
                  sig[2, 3] + (sig[1, 2]^2) * sig[2, 3] + (sig[1, 3]^2) * 
                  sig[2, 3] - (sig[2, 3]^3))))) / .Power(-1 + (sig[1, 2]^2) + 
                  (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * sig[2, 3] + 
                  (sig[2, 3]^2), 3)
  return(res.f)
}

.S.ga.12.13 <- function(u, sig) {
  res.f <- (-((-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * 
                 sig[2, 3] + (sig[2, 3]^2)) * (2 * sig[1, 2] * sig[1, 3] + 2 * 
                  (u[3]^2) * sig[1, 2] * sig[1, 3] + sig[2, 3] - (u[3]^2) * 
                  sig[2, 3] - 3 * (sig[1, 2]^2) * sig[2, 3] - (u[3]^2) * 
                  (sig[1, 2]^2) * sig[2, 3] - 3 * (sig[1, 3]^2) * sig[2, 3] + 
                  4 * sig[1, 2] * sig[1, 3] * (sig[2, 3]^2) - (sig[2, 3]^3) - 
                  2 * u[1] * u[3] * (sig[1, 2] - sig[1, 3] * sig[2, 3]) - 
                  (u[2]^2) * (2 * sig[1, 2] * sig[1, 3] + sig[2, 3] - 3 * 
                  (sig[1, 3]^2) * sig[2, 3]) + (u[1]^2) * sig[2, 3] * (-1 + 
                  (sig[2, 3]^2)) + u[2] * (u[3] * (1 + (sig[1, 2]^2) - 3 * 
                  (sig[1, 3]^2) + (sig[2, 3]^2)) + 2 * u[1] * (sig[1, 3] + 
                  sig[1, 2] * sig[2, 3] - 2 * sig[1, 3] * (sig[2, 3]^2))))) + 
                  4 * (sig[1, 3] - sig[1, 2] * sig[2, 3]) * ((sig[1, 2]^3) + 
                  (sig[1, 2]^2) * sig[1, 3] * (u[2] * u[3] - (3 + (u[3]^2)) * 
                  sig[2, 3]) - (u[1]^2) * (sig[1, 2] - sig[1, 3] * sig[2, 3]) * 
                  (-1 + (sig[2, 3]^2)) + sig[1, 3] * ((u[2]^2) * (-1 + 
                  (sig[1, 3]^2)) * sig[2, 3] + u[2] * u[3] * (1 - 
                  (sig[1, 3]^2) + (sig[2, 3]^2)) - sig[2, 3] * (-1 + (u[3]^2) + 
                  (sig[1, 3]^2) + (sig[2, 3]^2))) + sig[1, 2] * (-1 - 
                  (u[2]^2) * (-1 + (sig[1, 3]^2)) - 2 * u[2] * u[3] * 
                  sig[2, 3] + (sig[2, 3]^2) + (u[3]^2) * (sig[2, 3]^2) + 
                  (sig[1, 3]^2) * (1 + (u[3]^2) + 2 * (sig[2, 3]^2))) + 
                  u[1] * (u[2] * (-1 - (sig[1, 2]^2) + (sig[1, 3]^2) + 2 * 
                  sig[1, 2] * sig[1, 3] * sig[2, 3] + (sig[2, 3]^2) - 2 * 
                  (sig[1, 3]^2) * (sig[2, 3]^2)) + u[3] * (-2 * sig[1, 2] * 
                  sig[1, 3] + sig[2, 3] + (sig[1, 2]^2) * sig[2, 3] + 
                  (sig[1, 3]^2) * sig[2, 3] - (sig[2, 3]^3))))) / 
                  .Power(-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * 
                  sig[1, 3] * sig[2, 3] + (sig[2, 3]^2), 3)
  return(res.f)
}

.S.ga.12.23 <- function(u, sig) {
  res.f <- (-((-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * 
                 sig[1, 3] * sig[2, 3] + (sig[2, 3]^2)) * (-((3 + (u[3]^2)) * 
                  (sig[1, 2]^2) * sig[1, 3]) + 2 * (u[1]^2) * sig[2, 3] * 
                  (-sig[1, 2] + sig[1, 3] * sig[2, 3]) + 2 * sig[1, 2] * 
                  (-(u[2] * u[3]) + (1 + (u[3]^2) + 2 * (sig[1, 3]^2)) * 
                  sig[2, 3]) + sig[1, 3] * (1 - (u[3]^2) - (sig[1, 3]^2) + 
                  (u[2]^2) * (-1 + (sig[1, 3]^2)) + 2 * u[2] * u[3] * 
                  sig[2, 3] - 3 * (sig[2, 3]^2)) + (u[1]^2) * sig[1, 3] * 
                  (-1 + (sig[2, 3]^2)) + u[1] * (2 * u[2] * (sig[1, 2] * 
                  sig[1, 3] + sig[2, 3] - 2 * (sig[1, 3]^2) * sig[2, 3]) + 
                  u[3] * (1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 3 * 
                  (sig[2, 3]^2))))) + 2 * (-2 * sig[1, 2] * sig[1, 3] + 2 * 
                  sig[2, 3]) * ((sig[1, 2]^3) + (sig[1, 2]^2) * sig[1, 3] * 
                  (u[2] * u[3] - (3 + (u[3]^2)) * sig[2, 3]) - (u[1]^2) * 
                  (sig[1, 2] - sig[1, 3] * sig[2, 3]) * (-1 + (sig[2, 3]^2)) + 
                  sig[1, 3] * ((u[2]^2) * (-1 + (sig[1, 3]^2)) * sig[2, 3] + 
                  u[2] * u[3] * (1 - (sig[1, 3]^2) + (sig[2, 3]^2)) - 
                  sig[2, 3] * (-1 + (u[3]^2) + (sig[1, 3]^2) + 
                  (sig[2, 3]^2))) + sig[1, 2] * (-1 - (u[2]^2) * (-1 + 
                  (sig[1, 3]^2)) - 2 * u[2] * u[3] * sig[2, 3] + 
                  (sig[2, 3]^2) + (u[3]^2) * (sig[2, 3]^2) + (sig[1, 3]^2) * 
                  (1 + (u[3]^2) + 2 * (sig[2, 3]^2))) + u[1] * (u[2] * 
                  (-1 - (sig[1, 2]^2) + (sig[1, 3]^2) + 2 * sig[1, 2] * 
                  sig[1, 3] * sig[2, 3] + (sig[2, 3]^2) - 2 * (sig[1, 3]^2) * 
                  (sig[2, 3]^2)) + u[3] * (-2 * sig[1, 2] * sig[1, 3] + 
                  sig[2, 3] + (sig[1, 2]^2) * sig[2, 3] + (sig[1, 3]^2) * 
                  sig[2, 3] - (sig[2, 3]^3))))) / .Power(-1 + (sig[1, 2]^2) + 
                  (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * sig[2, 3] + 
                  (sig[2, 3]^2), 3)
  return(res.f)
}

.S.ga.13.13 <- function(u, sig) {
  res.f <- -(((-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * 
              sig[2, 3] + (sig[2, 3]^2)) * (-1 + (sig[1, 2]^2) + (u[2]^2) * 
              (sig[1, 2]^2) - (u[3]^2) * (-1 + (sig[1, 2]^2)) + 3 * 
              (sig[1, 3]^2) + 2 * u[2] * u[3] * (sig[1, 2] * sig[1, 3] - 
              sig[2, 3]) - 6 * sig[1, 2] * sig[1, 3] * sig[2, 3] - 2 * 
              (u[2]^2) * sig[1, 2] * sig[1, 3] * sig[2, 3] + (sig[2, 3]^2) + 
              (u[2]^2) * (sig[2, 3]^2) + 2 * (sig[1, 2]^2) * (sig[2, 3]^2) - 
              (u[1]^2) * (-1 + (sig[2, 3]^2)) - 2 * u[1] * (u[3] * (sig[1, 3] - 
              sig[1, 2] * sig[2, 3]) + u[2] * (sig[1, 2] - sig[1, 3] * 
              sig[2, 3]))) + 4 * (sig[1, 3] - sig[1, 2] * sig[2, 3]) * 
              ((u[1]^2) * (sig[1, 3] - sig[1, 2] * sig[2, 3]) * (-1 + 
              (sig[2, 3]^2)) - (sig[1, 3] - sig[1, 2] * sig[2, 3]) * (-1 + 
              (sig[1, 2]^2) - (u[3]^2) * (-1 + (sig[1, 2]^2)) + 
              (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * sig[2, 3] + 
              (sig[2, 3]^2)) - (u[2]^2) * ((sig[1, 2]^2) * sig[1, 3] - 
              sig[1, 2] * (1 + (sig[1, 3]^2)) * sig[2, 3] + sig[1, 3] * 
              (sig[2, 3]^2)) - u[2] * u[3] * (-(sig[1, 2]^3) - 2 * sig[1, 3] * 
              sig[2, 3] + sig[1, 2] * (1 + (sig[1, 3]^2) + (sig[2, 3]^2))) + 
              u[1] * (-(u[2] * (-2 * sig[1, 2] * sig[1, 3] + sig[2, 3] + 
              (sig[1, 2]^2) * sig[2, 3] + (sig[1, 3]^2) * sig[2, 3] - 
              (sig[2, 3]^3))) + u[3] * (1 + (sig[1, 3]^2) - 2 * sig[1, 2] * 
              sig[1, 3] * sig[2, 3] - (sig[2, 3]^2) + (sig[1, 2]^2) * 
              (-1 + 2 * (sig[2, 3]^2)))))) / .Power(-1 + (sig[1, 2]^2) + 
              (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * sig[2, 3] + 
              (sig[2, 3]^2), 3))
  return(res.f)
}

.S.ga.13.23 <- function(u, sig) {
  res.f <- ((-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * 
              sig[2, 3] + (sig[2, 3]^2)) * (-((-1 + (u[3]^2)) * 
              (sig[1, 2]^3)) - 4 * (sig[1, 2]^2) * sig[1, 3] * sig[2, 3] - 2 * 
              sig[1, 3] * (-(u[2] * u[3]) + sig[2, 3] + (u[2]^2) * sig[2, 3]) + 
              sig[1, 2] * (-1 + (u[3]^2) + 3 * (sig[1, 3]^2) + (u[2]^2) * 
              (1 + (sig[1, 3]^2)) - 2 * u[2] * u[3] * sig[2, 3] + 3 * 
              (sig[2, 3]^2)) + (u[1]^2) * (sig[1, 2] + 2 * sig[1, 3] * 
              sig[2, 3] - 3 * sig[1, 2] * (sig[2, 3]^2)) - u[1] * (2 * u[3] * 
              (sig[1, 2] * sig[1, 3] + sig[2, 3] - 2 * (sig[1, 2]^2) * 
              sig[2, 3]) + u[2] * (1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 3 * 
              (sig[2, 3]^2)))) - 2 * (-2 * sig[1, 2] * sig[1, 3] + 2 * 
              sig[2, 3]) * ((u[1]^2) * (sig[1, 3] - sig[1, 2] * sig[2, 3]) * 
              (-1 + (sig[2, 3]^2)) - (sig[1, 3] - sig[1, 2] * sig[2, 3]) * 
              (-1 + (sig[1, 2]^2) - (u[3]^2) * (-1 + (sig[1, 2]^2)) + 
              (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * sig[2, 3] + 
              (sig[2, 3]^2)) - (u[2]^2) * ((sig[1, 2]^2) * sig[1, 3] - 
              sig[1, 2] * (1 + (sig[1, 3]^2)) * sig[2, 3] + sig[1, 3] * 
              (sig[2, 3]^2)) - u[2] * u[3] * (-(sig[1, 2]^3) - 2 * 
              sig[1, 3] * sig[2, 3] + sig[1, 2] * (1 + (sig[1, 3]^2) + 
              (sig[2, 3]^2))) + u[1] * (-(u[2] * (-2 * sig[1, 2] * sig[1, 3] + 
              sig[2, 3] + (sig[1, 2]^2) * sig[2, 3] + (sig[1, 3]^2) * 
              sig[2, 3] - (sig[2, 3]^3))) + u[3] * (1 + (sig[1, 3]^2) - 2 * 
              sig[1, 2] * sig[1, 3] * sig[2, 3] - (sig[2, 3]^2) + 
              (sig[1, 2]^2) * (-1 + 2 * (sig[2, 3]^2)))))) / .Power(-1 + 
              (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * 
              sig[2, 3] + (sig[2, 3]^2), 3)
  return(res.f)
}

.S.ga.23.23 <- function(u, sig) {
  res.f <- (-((-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * sig[1, 3] * 
                sig[2, 3] + (sig[2, 3]^2)) * (-1 + (sig[1, 2]^2) + (u[1]^2) * 
                (sig[1, 2]^2) - (u[3]^2) * (-1 + (sig[1, 2]^2)) + 
                (sig[1, 3]^2) + (u[1]^2) * (sig[1, 3]^2) + 2 * (sig[1, 2]^2) * 
                (sig[1, 3]^2) - (u[2]^2) * (-1 + (sig[1, 3]^2)) - 6 * 
                sig[1, 2] * sig[1, 3] * sig[2, 3] - 2 * (u[1]^2) * sig[1, 2] * 
                sig[1, 3] * sig[2, 3] + 3 * (sig[2, 3]^2) - 2 * u[1] * u[3] * 
                (sig[1, 3] - sig[1, 2] * sig[2, 3]) + 2 * u[2] * (-(u[1] * 
                sig[1, 2]) + u[3] * sig[1, 2] * sig[1, 3] - u[3] * sig[2, 3] + 
                u[1] * sig[1, 3] * sig[2, 3]))) + 2 * (-2 * sig[1, 2] * 
                sig[1, 3] + 2 * sig[2, 3]) * ((u[2]^2) * (-1 + (sig[1, 3]^2)) * 
                (sig[1, 2] * sig[1, 3] - sig[2, 3]) + (sig[1, 2] * sig[1, 3] - 
                sig[2, 3]) * (1 - (sig[1, 2]^2) + (u[3]^2) * 
                (-1 + (sig[1, 2]^2)) - (sig[1, 3]^2) + 2 * sig[1, 2] * 
                sig[1, 3] * sig[2, 3] - (sig[2, 3]^2)) + (u[1]^2) * 
                ((sig[1, 2]^2) * sig[2, 3] + (sig[1, 3]^2) * sig[2, 3] - 
                sig[1, 2] * sig[1, 3] * (1 + (sig[2, 3]^2))) + u[1] * u[3] * 
                (-(sig[1, 2]^3) - 2 * sig[1, 3] * sig[2, 3] + sig[1, 2] * 
                (1 + (sig[1, 3]^2) + (sig[2, 3]^2))) + u[2] * (u[3] * (-1 + 
                (sig[1, 3]^2) + (sig[1, 2]^2) * (1 - 2 * (sig[1, 3]^2)) + 2 * 
                sig[1, 2] * sig[1, 3] * sig[2, 3] - (sig[2, 3]^2)) + u[1] * 
                (-(sig[1, 3]^3) - 2 * sig[1, 2] * sig[2, 3] + sig[1, 3] * 
                (1 + (sig[1, 2]^2) + (sig[2, 3]^2)))))) / 
                .Power(-1 + (sig[1, 2]^2) + (sig[1, 3]^2) - 2 * sig[1, 2] * 
                sig[1, 3] * sig[2, 3] + (sig[2, 3]^2), 3)
  return(res.f)
}