R/CDVineSim.R

In elong0527/simuCPP: Simulate Concordance Probability, AUC and Time-Dependent AUC Based on Copula

#' Function from CDVine
#' https://github.com/cran/CDVine
#' @noRd
CDVineSim <- function(N, family, par, par2 = rep(0, length(family)), type) {
  ## #############################################################################
  ##  Function that simulates copulae, gaussian, student's t, clayton or gumbel ##
  ## -------------------------------------------------------------------------- ##
  ##  INPUT:                                                                    ##
  ##    N                            Number of simulations                      ##
  ##    d                            Dimension of copula                        ##
  ##    family                       Array defining the ppc                     ##
  ##    par  		       Dependency parameters                      ##
  ##    par2			       Second copula parameter			  ##
  ##	  type			       Vine type (1=C-vine, 2=D-vine)		  ##
  ##  OUTPUT:                                                                   ##
  ##    U                            Simulated dxN copula                       ##
  ## -------------------------------------------------------------------------- ##
  ## #############################################################################

  if (type == "CVine")
    type <- 1 else if (type == "DVine")
      type <- 2
    if (type != 1 & type != 2)
      stop("Vine model not implemented.")

    dd <- length(family)
    d <- (1 + sqrt(1 + 8 * dd))/2
    if (d != floor(d))
      stop("The length of the family vector is not correct.")

    if (length(par) < dd)
      stop("Length of 'par' has to be d*(d-1)/2.")
    if (any(c(2, 7, 8, 9, 10, 17, 18, 19, 20, 27, 28, 29, 30, 37, 38, 39, 40) %in% family) && length(par2) != dd)
      stop("Length of 'par2' has to be d*(d-1)/2 if t, BB1, BB6, BB7 or BB8 copulas are used.")
    if (length(family) != dd)
      stop("Length of 'family' has to be d*(d-1)/2.")
    if (d < 2)
      stop("Dimension has to be at least 2.")

    for (i in 1:(d * (d - 1)/2)) {
      if (!(family[i] %in% c(0, 1:10, 13, 14, 16:20, 23, 24, 26:30, 33, 34, 36:40)))
        stop("Copula family not implemented.")
      # Parameterbereiche abfragen
      if ((family[i] == 1 || family[i] == 2) && abs(par[i]) >= 1)
        stop("The parameter of the Gaussian and t-copula has to be in the interval (-1,1).")
      if (family[i] == 2 && par2[i] <= 2)
        stop("The degrees of freedom parameter of the t-copula has to be larger than 2.")
      if ((family[i] == 3 || family[i] == 13) && par[i] <= 0)
        stop("The parameter of the Clayton copula has to be positive.")
      if ((family[i] == 4 || family[i] == 14) && par[i] < 1)
        stop("The parameter of the Gumbel copula has to be in the interval [1,oo).")
      if ((family[i] == 6 || family[i] == 16) && par[i] <= 1)
        stop("The copula parameter of the Joe copula has to be in the interval (1,oo).")
      if (family[i] == 5 && par[i] == 0)
        stop("The parameter of the Frank copula has to be unequal to 0.")
      if ((family[i] == 7 || family[i] == 17) && par[i] <= 0)
        stop("The first parameter of the BB1 copula has to be positive.")
      if ((family[i] == 7 || family[i] == 17) && par2[i] < 1)
        stop("The second parameter of the BB1 copula has to be in the interval [1,oo).")
      if ((family[i] == 8 || family[i] == 18) && par[i] <= 0)
        stop("The first parameter of the BB6 copula has to be in the interval [1,oo).")
      if ((family[i] == 8 || family[i] == 18) && par2[i] < 1)
        stop("The second parameter of the BB6 copula has to be in the interval [1,oo).")
      if ((family[i] == 9 || family[i] == 19) && par[i] < 1)
        stop("The first parameter of the BB7 copula has to be in the interval [1,oo).")
      if ((family[i] == 9 || family[i] == 19) && par2[i] <= 0)
        stop("The second parameter of the BB7 copula has to be positive.")
      if ((family[i] == 10 || family[i] == 20) && par[i] < 1)
        stop("The first parameter of the BB8 copula has to be in the interval [1,oo).")
      if ((family[i] == 10 || family[i] == 20) && (par2[i] <= 0 || par2[i] > 1))
        stop("The second parameter of the BB8 copula has to be in the interval (0,1].")
      if ((family[i] == 23 || family[i] == 33) && par[i] >= 0)
        stop("The parameter of the rotated Clayton copula has to be negative.")
      if ((family[i] == 24 || family[i] == 34) && par[i] > -1)
        stop("The parameter of the rotated Gumbel copula has to be in the interval (-oo,-1].")
      if ((family[i] == 26 || family[i] == 36) && par[i] >= -1)
        stop("The parameter of the rotated Joe copula has to be in the interval (-oo,-1).")
      if ((family[i] == 27 || family[i] == 37) && par[i] >= 0)
        stop("The first parameter of the rotated BB1 copula has to be negative.")
      if ((family[i] == 27 || family[i] == 37) && par2[i] > -1)
        stop("The second parameter of the rotated BB1 copula has to be in the interval (-oo,-1].")
      if ((family[i] == 28 || family[i] == 38) && par[i] >= 0)
        stop("The first parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
      if ((family[i] == 28 || family[i] == 38) && par2[i] > -1)
        stop("The second parameter of the rotated BB6 copula has to be in the interval (-oo,-1].")
      if ((family[i] == 29 || family[i] == 39) && par[i] > -1)
        stop("The first parameter of the rotated BB7 copula has to be in the interval (-oo,-1].")
      if ((family[i] == 29 || family[i] == 39) && par2[i] >= 0)
        stop("The second parameter of the rotated BB7 copula has to be negative.")
      if ((family[i] == 30 || family[i] == 40) && par[i] > -1)
        stop("The first parameter of the rotated BB8 copula has to be in the interval (-oo,-1].")
      if ((family[i] == 30 || family[i] == 40) && (par2[i] >= 0 || par2[i] < (-1)))
        stop("The second parameter of the rotated BB8 copula has to be in the interval [-1,0).")
    }

    if (type == 1) {
      if (!any(family %in% c(2, 7:10, 17:20, 27:30, 37:40)))
        tmp <- .C("pcc", as.integer(N), as.integer(d), as.integer(family), as.integer(type), as.double(par),
                  as.double(rep(0, d * (d - 1)/2)), as.double(rep(0, N * d)), PACKAGE = "simuCCP")[[7]] else tmp <- .C("pcc", as.integer(N), as.integer(d), as.integer(family), as.integer(type), as.double(par),
                                                                                                                      as.double(par2), as.double(rep(0, N * d)), PACKAGE = "simuCCP")[[7]]
                  U <- matrix(tmp, ncol = d)

    } else if (type == 2) {
      # D-vine

      Matrix <- matrix(rep(0, d * d), d, d)
      Copula.Params <- matrix(rep(0, d * d), d, d)
      Copula.Params2 <- matrix(rep(0, d * d), d, d)
      Copula.Types <- matrix(rep(0, d * d), d, d)

      for (i in 1:d) {
        Matrix[(d - i + 1), (d - i + 1)] <- i
      }

      k <- 1
      for (i in 1:(d - 1)) {
        for (j in 1:(d - i)) {
          Matrix[(d - i + 1), (d - j - i + 1)] <- j
          Copula.Types[(d - i + 1), (d - j - i + 1)] <- family[k]
          Copula.Params[(d - i + 1), (d - j - i + 1)] <- par[k]
          Copula.Params2[(d - i + 1), (d - j - i + 1)] <- par2[k]
          k <- k + 1
        }
      }

      Matrix[is.na(Matrix)] <- 0
      Copula.Types[is.na(Copula.Types)] <- 0
      Copula.Params[is.na(Copula.Params)] <- 0
      Copula.Params2[is.na(Copula.Params2)] <- 0

      MaxMat <- createMaxMat(Matrix)
      CondDistr <- neededCondDistr(Matrix)

      dvlist <- list(Matrix = Matrix, family = Copula.Types, par = Copula.Params, par2 = Copula.Params2,
                     MaxMat = MaxMat, CondDistr = CondDistr)

      o <- diag(dvlist$Matrix)
      dvlist <- normalizeDv(dvlist)

      matri <- as.vector(dvlist$Matrix)
      w1 <- as.vector(dvlist$family)
      th <- as.vector(dvlist$par)
      th2 <- as.vector(dvlist$par2)
      maxmat <- as.vector(dvlist$MaxMat)
      conindirect <- as.vector(dvlist$CondDistr$indirect)
      matri[is.na(matri)] <- 0
      w1[is.na(w1)] <- 0
      th[is.na(th)] <- 0
      th2[is.na(th2)] <- 0
      maxmat[is.na(maxmat)] <- 0
      conindirect[is.na(conindirect)] <- 0

      tmp <- rep(0, d * N)

      tmp <- .C("SimulateVine", as.integer(N), as.integer(d), as.integer(w1), as.integer(maxmat), as.integer(matri),
                as.integer(conindirect), as.double(th), as.double(th2), as.double(tmp), PACKAGE = "simuCCP")[[9]]

      out <- matrix(tmp, ncol = d)


      U <- out[, sort(o[length(o):1], index.return = TRUE)$ix]

    }
    return(U)
}

normalizeDv <- function(dvlist) {

  oldOrder <- diag(dvlist$Matrix)
  Matrix <- reorderDv(dvlist$Matrix)

  dvlist$Matrix[is.na(dvlist$Matrix)] <- 0
  dvlist$family[is.na(dvlist$family)] <- 0
  dvlist$par[is.na(dvlist$par)] <- 0
  dvlist$par2[is.na(dvlist$par2)] <- 0

  MaxMat <- createMaxMat(Matrix)
  CondDistr <- neededCondDistr(Matrix)

  return(list(Matrix = Matrix, family = dvlist$family, par = dvlist$par, par2 = dvlist$par2, MaxMat = MaxMat,
              CondDistr = CondDistr))
}

reorderDv <- function(Matrix) {
  oldOrder <- diag(Matrix)

  O <- apply(t(1:nrow(Matrix)), 2, "==", Matrix)

  for (i in 1:nrow(Matrix)) {
    Matrix[O[, oldOrder[i]]] <- nrow(Matrix) - i + 1
  }

  return(Matrix)
}

createMaxMat <- function(Matrix) {

  if (dim(Matrix)[1] != dim(Matrix)[2])
    stop("Structure matrix has to be quadratic.")

  MaxMat <- reorderDv(Matrix)

  n <- nrow(MaxMat)

  for (j in 1:(n - 1)) {
    for (i in (n - 1):j) {
      MaxMat[i, j] <- max(MaxMat[i:(i + 1), j])
    }
  }

  tMaxMat <- MaxMat
  tMaxMat[is.na(tMaxMat)] <- 0

  oldSort <- diag(Matrix)
  oldSort <- oldSort[n:1]

  for (i in 1:n) {
    MaxMat[tMaxMat == i] <- oldSort[i]
  }

  return(MaxMat)
}

neededCondDistr <- function(Vine) {

  if (dim(Vine)[1] != dim(Vine)[2])
    stop("Structure matrix has to be quadratic.")

  Vine <- reorderDv(Vine)

  MaxMat <- createMaxMat(Vine)

  d <- nrow(Vine)

  M <- list()
  M$direct <- matrix(FALSE, d, d)
  M$indirect <- matrix(FALSE, d, d)

  M$direct[2:d, 1] <- TRUE

  for (i in 2:(d - 1)) {
    v <- d - i + 1

    bw <- as.matrix(MaxMat[i:d, 1:(i - 1)]) == v

    direct <- Vine[i:d, 1:(i - 1)] == v

    M$indirect[i:d, i] <- apply(as.matrix(bw & (!direct)), 1, any)

    M$direct[i:d, i] <- TRUE

    M$direct[i, i] <- any(as.matrix(bw)[1, ] & as.matrix(direct)[1, ])
  }

  return(M)
}