R/pExtDep.R

In ExtremalDep: Extremal Dependence Models

pExtDep <- function(q, type, method = "Parametric", model, par, plot = TRUE,
                    main, xlab, cex.lab, cex.axis, lwd, ...) {
  # Checking quantile vector or matrix
  if (is.vector(q)) {
    dim <- length(q)
    q0 <- q
    if (is.matrix(par)) {
      npar <- nrow(par)
      p <- vector(length = npar)
      if (any(q == 0)) {
        return(rep(NA, npar))
      }
    } else if (is.vector(par)) {
      if (any(q == 0)) {
        return(NA)
      }
    }
  } else if (is.matrix(q)) {
    dim <- ncol(q)
    qNA <- apply(q, 1, function(x) any(x == 0))

    q0 <- q
    nq0 <- nrow(q0)
    q <- q[!qNA, ]
    if (is.matrix(q)) {
      nq <- nrow(q)
    } else if (is.vector(q)) {
      nq <- 1
    }

    if (is.vector(par)) {
      p <- vector(length = nq)
    } else if (is.matrix(par)) {
      p <- matrix(nrow = nrow(par), ncol = nq)
    }
  } else {
    stop("q must be a vector or a matrix")
  }

  # Checking method
  methods <- c("Parametric", "NonParametric")
  if (!any(method == methods)) {
    stop("Wrong method specified")
  }

  # Checking the type
  types <- c("lower", "inv.lower", "upper")
  if (method == "Parametric" && !any(type == types)) {
    stop("Wrong probability type specified")
  }

  if (method == "Parametric") {
    if (all(dim != c(2, 3))) {
      stop("Dimension 2 or 3 only")
    }

    # Checking the model
    models <- c("HR", "ET", "EST")
    if (!any(model == models)) {
      stop("model wrongly specified")
    }

    # Checking model, par & dim are correctly specified
    if (is.vector(par)) {
      if (!dim_ExtDep(model = model, par = par, dim = dim)) {
        stop("Length of 'par' is incorrect")
      }
    } else if (is.matrix(par)) {
      if (!dim_ExtDep(model = model, par = par[1, ], dim = dim)) {
        stop("Length of 'par' is incorrect")
      }
    } else {
      stop("'par' should be a vector or a matrix")
    }

    if (model == "HR") {
      if (is.vector(q)) {
        if (is.vector(par)) {
          p <- p.hr(q = q, par = par, type = type)
        } else if (is.matrix(par)) {
          p <- apply(par, 1, p.hr, q = q, type = type)
        }
      } else if (is.matrix(q)) {
        if (is.vector(par)) {
          p <- apply(q, 1, p.hr, par = par, type = type)
        } else if (is.matrix(par)) {
          for (i in 1:nq) {
            p[, i] <- apply(par, 1, p.hr, q = q[i, ], type = type)
          }
        }
      }
    } else if (model == "ET") {
      if (is.vector(q)) {
        if (is.vector(par)) {
          p <- p.et(q = q, par = par, type = type)
        } else if (is.matrix(par)) {
          p <- apply(par, 1, p.et, q = q, type = type)
        }
      } else if (is.matrix(q)) {
        if (is.vector(par)) {
          p <- apply(q, 1, p.et, par = par, type = type)
        } else if (is.matrix(par)) {
          for (i in 1:nq) {
            p[, i] <- apply(par, 1, p.et, q = q[i, ], type = type)
          }
        }
      }
    } else if (model == "EST") {
      if (is.vector(q)) {
        if (is.vector(par)) {
          p <- p.est(q = q, par = par, type = type)
        } else if (is.matrix(par)) {
          p <- apply(par, 1, p.est, q = q, type = type)
        }
      } else if (is.matrix(q)) {
        if (is.vector(par)) {
          p <- apply(q, 1, p.est, par = par, type = type)
        } else if (is.matrix(par)) {
          for (i in 1:nq) {
            p[, i] <- apply(par, 1, p.est, q = q[i, ], type = type)
          }
        }
      }
    }
  } else if (method == "NonParametric") {
    if (is.vector(q)) {
      if (is.vector(par)) {
        p <- ph(w = q, beta = par)
      } else if (is.matrix(par)) {
        for (i in 1:npar) {
          p[i] <- ph(w = q, beta = par[i, ])
        }
      }
    } else if (is.matrix(q)) {
      if (is.vector(par)) {
        for (i in 1:nq) {
          p[i] <- ph(w = q[i, ], beta = par)
        }
      } else if (is.matrix(par)) {
        for (i in 1:nq) {
          for (j in 1:npar) {
            p[j, i] <- ph(w = q[i, ], beta = par[j, ])
          }
        }
      }
    }
  }

  if (is.matrix(par)) {
    if (plot) {
      if (missing(main)) {
        main <- ""
      }
      if (missing(cex.lab)) {
        cex.lab <- 1.4
      }
      if (missing(cex.axis)) {
        cex.axis <- 1.4
      }
      if (missing(lwd)) {
        lwd <- 2
      }

      if (is.vector(p)) {
        if (missing(xlab)) {
          if (dim == 2) {
            if (type == "lower") {
              xlab <- paste("P(X<", q[1], ",Y<", q[2], ")", sep = "")
            } else if (type == "inv.lower") {
              xlab <- paste("1-P(X<", q[1], ",Y<", q[2], ")", sep = "")
            } else if (type == "upper") {
              xlab <- paste("P(X>", q[1], ",Y>", q[2], ")", sep = "")
            }
          } else if (dim == 3) {
            if (type == "lower") {
              xlab <- paste("P(X<", q[1], ",Y<", q[2], ",Z<", q[3], ")", sep = "")
            } else if (type == "inv.lower") {
              xlab <- paste("1-P(X<", q[1], ",Y<", q[2], ",Z<", q[3], ")", sep = "")
            } else if (type == "upper") {
              xlab <- paste("P(X>", q[1], ",Y>", q[2], ",Z>", q[3], ")", sep = "")
            }
          }
        }

        Ke <- density(p) # KDE
        Hi <- hist(p,
          prob = TRUE, col = "lightgrey",
          ylim = range(Ke$y), main = main, xlab = xlab,
          cex.lab = cex.lab, cex.axis = cex.axis, lwd = lwd, ...
        )
        p_ic <- quantile(p, probs = c(0.025, 0.5, 0.975))

        points(x = p_ic, y = c(0, 0, 0), pch = 4, lwd = 4)
        points(x = mean(p), y = 0, pch = 16, lwd = 4)
        lines(Ke, lwd = 2, col = "dimgrey")
      } else if (is.matrix(p)) {
        for (i in 1:nq) {
          if (missing(xlab)) {
            if (dim == 2) {
              if (type == "lower") {
                xlab <- paste("P(X<", q[i, 1], ",Y<", q[i, 2], ")", sep = "")
              } else if (type == "inv.lower") {
                xlab <- paste("1-P(X<", q[i, 1], ",Y<", q[i, 2], ")", sep = "")
              } else if (type == "upper") {
                xlab <- paste("P(X>", q[i, 1], ",Y>", q[i, 2], ")", sep = "")
              }
            } else if (dim == 3) {
              if (type == "lower") {
                xlab <- paste("P(X<", q[i, 1], ",Y<", q[i, 2], ",Z<", q[i, 3], ")", sep = "")
              } else if (type == "inv.lower") {
                xlab <- paste("1-P(X<", q[i, 1], ",Y<", q[i, 2], ",Z<", q[i, 3], ")", sep = "")
              } else if (type == "upper") {
                xlab <- paste("P(X>", q[i, 1], ",Y>", q[i, 2], ",Z>", q[i, 3], ")", sep = "")
              }
            }
          }

          Ke <- density(p[, i]) # KDE
          Hi <- hist(p[, i],
            prob = TRUE, col = "lightgrey",
            ylim = range(Ke$y), main = main, xlab = xlab,
            cex.lab = cex.lab, cex.axis = cex.axis, lwd = lwd, ...
          )
          p_ic <- quantile(p[, i], probs = c(0.025, 0.5, 0.975))

          points(x = p_ic, y = c(0, 0, 0), pch = 4, lwd = 4)
          points(x = mean(p[, i]), y = 0, pch = 16, lwd = 4)
          lines(Ke, lwd = 2, col = "dimgrey")
        }
      }
    }
  }

  if (is.matrix(q0)) {
    if (sum(qNA) != 0) {
      p.temp <- p

      if (is.vector(par)) {
        p <- vector(length = nq0)
        p[!qNA] <- p.temp
        p[qNA] <- rep(NA, sum(qNA))
      } else if (is.matrix(par)) {
        p <- matrix(nrow = nrow(par), ncol = nq0)
        p[, !qNA] <- p.temp
        # p[,qNA] <- matrix(NA, nrow=nrow(par), ncol=sum(qNA))
      }
    }
  }

  return(p)
}

pFailure <- function(n, beta, u1, u2, mar1, mar2, type, plot, xlab, ylab, nlevels = 10) {
  if (!(type %in% c("and", "or", "both"))) {
    stop(" 'type' must be 'and', 'or' or 'both'.")
  }

  ### Simulate data

  w <- rh(n, beta)
  r <- runif(n)^(-1)
  y <- cbind(2 * r * w, 2 * r * (1 - w))

  ### Compute ustar

  u.star <- function(threshold, mar) {
    if (mar[3] == 0) {
      return(exp((threshold - mar[1]) / mar[2]))
    } else {
      return((1 + mar[3] * (threshold - mar[1]) / mar[2])^(1 / mar[3]))
    }
  }

  ustar1 <- u.star(u1, mar1)
  ustar2 <- u.star(u2, mar2)

  ustar <- expand.grid(ustar1, ustar2)

  if (type == "and") {
    AND <- apply(ustar, 1, function(u) sum(y[, 1] > u[1] & y[, 2] > u[2])) / n
    AND <- matrix(AND, nrow = length(u1), ncol = length(u2), byrow = FALSE)
  }
  if (type == "or") {
    OR <- apply(ustar, 1, function(u) sum(y[, 1] > u[1] | y[, 2] > u[2])) / n
    OR <- matrix(OR, nrow = length(u1), ncol = length(u2), byrow = FALSE)
  }
  if (type == "both") {
    and_or <- function(u, y) {
      c(sum(y[, 1] > u[1] & y[, 2] > u[2]), sum(y[, 1] > u[1] | y[, 2] > u[2]))
    }

    AND_OR <- apply(ustar, 1, and_or, y = y) / n
    AND <- matrix(AND_OR[1, ], nrow = length(u1), ncol = length(u2), byrow = FALSE)
    OR <- matrix(AND_OR[2, ], nrow = length(u1), ncol = length(u2), byrow = FALSE)
  }


  if (plot) {
    op <- par(mai = c(0.8, 0.8, 0.35, 0.1), mgp = c(2.5, 1, 0), cex.axis = 2, cex.lab = 2, cex.main = 2)
    on.exit(par(op))

    if (type == "and") {
      contour(u1, u2, AND,
        levels = round(seq(min(AND), max(AND), length = nlevels), 3),
        xlab = xlab, ylab = ylab, main = "P(X>x and Y>y)", labcex = 1
      )
    }
    if (type == "or") {
      contour(u1, u2, OR,
        levels = round(seq(min(OR), max(OR), length = nlevels), 3),
        xlab = xlab, ylab = ylab, main = "P(X>x or Y>y)", labcex = 1
      )
    }
    if (type == "both") {
      contour(u1, u2, AND,
        levels = round(seq(min(AND), max(AND), length = nlevels), 3),
        xlab = xlab, ylab = ylab, main = "P(X>x and Y>y)", labcex = 1
      )

      contour(u1, u2, OR,
        levels = round(seq(min(OR), max(OR), length = nlevels), 3),
        xlab = xlab, ylab = ylab, main = "P(X>x or Y>y)", labcex = 1
      )
    }
  }

  if (type == "and") {
    return(list(AND = AND))
  }
  if (type == "or") {
    return(list(OR = OR))
  }
  if (type == "both") {
    return(list(AND = AND, OR = OR))
  }
}

#########################################
#########################################
### Internal functions for pExtDep
#########################################
#########################################

## Husler-Reiss model

p.hr <- function(q, par, type) {
  p.hr.2d <- function(q, par, type) {
    I1 <- pnorm(par + log(q[2] / q[1]) / (2 * par))
    I2 <- pnorm(par + log(q[1] / q[2]) / (2 * par))

    if (type == "lower") {
      return(exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    } else if (type == "inv.lower") {
      return(1 - exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    } else if (type == "upper") {
      return(1 - exp(-1 / q[1]) - exp(-1 / q[2]) + exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    }
  }

  p.hr.3d <- function(q, par, type) {
    lambda12 <- par[1]
    lambda13 <- par[2]
    lambda23 <- par[3]
    S1 <- matrix(c(4 * lambda12^2, rep(2 * (lambda12^2 + lambda13^2 - lambda23^2), 2), 4 * lambda13^2), nrow = 2)
    S2 <- matrix(c(4 * lambda12^2, rep(2 * (lambda12^2 + lambda23^2 - lambda13^2), 2), 4 * lambda23^2), nrow = 2)
    S3 <- matrix(c(4 * lambda13^2, rep(2 * (lambda13^2 + lambda23^2 - lambda12^2), 2), 4 * lambda23^2), nrow = 2)

    I1 <- pmesn(x = c(log(q[2] / q[1]) + 2 * lambda12^2, log(q[3] / q[1]) + 2 * lambda13^2), scale = S1)
    I2 <- pmesn(x = c(log(q[1] / q[2]) + 2 * lambda12^2, log(q[3] / q[2]) + 2 * lambda23^2), scale = S2)
    I3 <- pmesn(x = c(log(q[1] / q[3]) + 2 * lambda13^2, log(q[2] / q[3]) + 2 * lambda23^2), scale = S3)

    if (type == "lower") {
      return(exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    } else if (type == "inv.lower") {
      return(1 - exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    } else if (type == "upper") {
      return(1 - exp(-1 / q[1]) - exp(-1 / q[2]) - exp(-1 / q[3])
        + p.hr.2d(q[1:2], lambda12, type = "lower") + p.hr.2d(q[c(1, 3)], lambda13, type = "lower")
        + p.hr.2d(q[2:3], lambda23, type = "lower") - exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    }
  }

  if (any(par <= 0)) {
    stop("wrong value of parameters")
  }
  dim <- length(q)
  if (!dim_ExtDep(model = "HR", par = par, dim = dim)) {
    stop("Length of par is incorrect")
  }
  if (dim == 2) {
    return(p.hr.2d(q = q, par = par, type = type))
  } else if (dim == 3) {
    return(p.hr.3d(q = q, par = par, type = type))
  }
}

## Extremal-t model

p.et <- function(q, par, type) {
  p.et.2d <- function(q, par, type) {
    if (((par[1] <= -1) || (par[1] >= 1) || (par[2] <= 0)) == 1) {
      stop("Extremal-t parameters ill defined")
    }
    rho <- par[1]
    nu <- par[2]

    I1 <- pt(sqrt((nu + 1) / (1 - rho^2)) * ((q[2] / q[1])^(1 / nu) - rho), df = nu + 1)
    I2 <- pt(sqrt((nu + 1) / (1 - rho^2)) * ((q[1] / q[2])^(1 / nu) - rho), df = nu + 1)

    if (type == "lower") {
      return(exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    } else if (type == "inv.lower") {
      return(1 - exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    } else if (type == "upper") {
      return(1 / q[1] * (1 - I1) + 1 / q[2] * (1 - I2))
    }
  }

  p.et.3d <- function(q, par, type) {
    if ((any(par[1:3] <= -1) || any(par[1:3] >= 1) || (par[4] <= 0)) == 1) {
      stop("Extremal-t parameters ill defined")
    }

    rho12 <- par[1]
    rho13 <- par[2]
    rho23 <- par[3]
    nu <- par[4]
    R1 <- (rho23 - rho12 * rho13) / sqrt((1 - rho12^2) * (1 - rho13^2))
    R2 <- (rho13 - rho12 * rho23) / sqrt((1 - rho12^2) * (1 - rho23^2))
    R3 <- (rho12 - rho13 * rho23) / sqrt((1 - rho13^2) * (1 - rho23^2))
    if (sum(eigen(matrix(c(1, R1, R1, 1), ncol = 2))$values < 0) >= 1) {
      stop("negative definite partial correlation matrix")
    }
    if (sum(eigen(matrix(c(1, R2, R2, 1), ncol = 2))$values < 0) >= 1) {
      stop("negative definite partial correlation matrix")
    }
    if (sum(eigen(matrix(c(1, R3, R3, 1), ncol = 2))$values < 0) >= 1) {
      stop("negative definite partial correlation matrix")
    }

    x1 <- sqrt((nu + 1) / (1 - rho12^2)) * ((q[2] / q[1])^(1 / nu) - rho12)
    y1 <- sqrt((nu + 1) / (1 - rho13^2)) * ((q[3] / q[1])^(1 / nu) - rho13)
    x2 <- sqrt((nu + 1) / (1 - rho12^2)) * ((q[1] / q[2])^(1 / nu) - rho12)
    y2 <- sqrt((nu + 1) / (1 - rho23^2)) * ((q[3] / q[2])^(1 / nu) - rho23)
    x3 <- sqrt((nu + 1) / (1 - rho13^2)) * ((q[1] / q[3])^(1 / nu) - rho13)
    y3 <- sqrt((nu + 1) / (1 - rho23^2)) * ((q[2] / q[3])^(1 / nu) - rho23)

    I1 <- pmest(x = c(x1, y1), scale = matrix(c(1, R1, R1, 1), ncol = 2), df = nu + 1)
    I2 <- pmest(x = c(x2, y2), scale = matrix(c(1, R2, R2, 1), ncol = 2), df = nu + 1)
    I3 <- pmest(x = c(x3, y3), scale = matrix(c(1, R3, R3, 1), ncol = 2), df = nu + 1)

    if (type == "lower") {
      return(exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    } else if (type == "inv.lower") {
      return(1 - exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    } else if (type == "upper") {
      return(1 - exp(-1 / q[1]) - exp(-1 / q[2]) - exp(-1 / q[3])
        + p.et.2d(q[1:2], c(rho12, nu), type = "lower") + p.et.2d(q[c(1, 3)], c(rho13, nu), type = "lower")
        + p.et.2d(q[2:3], c(rho23, nu), type = "lower") - exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    }
  }

  dim <- length(q)
  if (!dim_ExtDep(model = "ET", par = par, dim = dim)) {
    stop("Length of par is incorrect")
  }
  if (dim == 2) {
    return(p.et.2d(q = q, par = par, type = type))
  } else if (dim == 3) {
    return(p.et.3d(q = q, par = par, type = type))
  }
}

## Extremal Skew-t

pmextst <- function(x, scale = 1, shape = rep(0, length(x)), df = 1) {
  if (any(is.na(x))) {
    return(NA)
  }
  .C("pmextst", as.double(x), as.double(scale), as.double(df),
    as.double(shape),
    out = double(1), NAOK = TRUE
  )$out
}

p.est <- function(q, par, type) {
  ## Preliminary functions

  Sigma <- function(rho) {
    p <- round(uniroot(function(x) {
      length(rho) - choose(x, 2)
    }, lower = 1, upper = 10)$root)
    Sig <- diag(p)
    Sig[lower.tri(Sig)] <- rho
    Sig[row(Sig) < col(Sig)] <- t(Sig)[row(Sig) < col(Sig)]
    return(Sig)
  }

  Sigma_j <- function(rho, j) {
    Sig <- Sigma(rho)
    return(Sig[-j, -j] - tcrossprod(Sig[-j, j], Sig[j, -j]))
  }

  s_j <- function(rho, j) {
    p <- round(uniroot(function(x) {
      length(rho) - choose(x, 2)
    }, lower = 1, upper = 10)$root)
    k <- p - 1

    sigma_j <- Sigma_j(rho, j)
    M <- diag(k)
    diag(M) <- sqrt(diag(sigma_j))
    return(M)
  }

  Sigma_bar_j <- function(rho, j) {
    sigma_j <- Sigma_j(rho, j)
    sj <- s_j(rho, j)
    sj_inv <- chol2inv(chol(sj))
    return(tcrossprod(sj_inv, sigma_j) %*% sj_inv)
  }

  alpha_tilde <- function(alpha, j) {
    return(t(alpha[-j]))
  }

  alpha_bar_j <- function(rho, alpha, j) {
    Sig <- Sigma(rho)
    sigma_j <- Sigma_j(rho, j)
    Alpha_tilde <- alpha_tilde(alpha, j)

    num <- alpha[j] + tcrossprod(Sig[j, -j], Alpha_tilde)
    denom <- sqrt(1 + tcrossprod(tcrossprod(Alpha_tilde, sigma_j), Alpha_tilde))
    return(num / denom)
  }

  alpha_star_j <- function(rho, alpha, j) {
    Alpha_tilde <- alpha_tilde(alpha, j)
    sj <- s_j(rho, j)
    return(Alpha_tilde %*% sj)
  }

  tau_star_j <- function(rho, alpha, nu, j) {
    Sig <- Sigma(rho)
    Alpha_tilde <- alpha_tilde(alpha, j)
    return(sqrt(nu + 1) * (alpha[j] + tcrossprod(Sig[-j, j], Alpha_tilde)))
  }

  nu_p <- function(rho, alpha, nu, j) {
    Alpha_bar <- alpha_bar_j(rho, alpha, j)
    return(pt(sqrt(nu + 1) * Alpha_bar, df = nu + 1) * 2^(nu / 2 - 1) * gamma((nu + 1) / 2) / sqrt(pi))
  }

  x_bar <- function(x, rho, alpha, nu, j) {
    return(x * nu_p(rho, alpha, nu, j))
  }


  p.est.2d <- function(q, par, type) {
    dim <- length(q)
    if (dim != 2) {
      stop("Dimension should be 2")
    }

    rho <- par[1]
    alpha <- par[2:3]
    nu <- par[4]

    if (((rho <= -1) || (rho >= 1) || (nu <= 0)) == 1) {
      stop("Extremal-t parameters ill defined")
    }

    nu_1 <- nu_p(rho = rho, alpha = alpha, nu = nu, j = 1)
    nu_2 <- nu_p(rho = rho, alpha = alpha, nu = nu, j = 2)

    x1_bar <- q[1] * nu_1
    x2_bar <- q[2] * nu_2

    comp1_1 <- sqrt((nu + 1) / (1 - rho^2)) * ((x2_bar / x1_bar)^(1 / nu) - rho)
    comp1_2 <- sqrt((nu + 1) / (1 - rho^2)) * ((x1_bar / x2_bar)^(1 / nu) - rho)

    alpha_star_1 <- alpha_star_j(rho = rho, alpha = alpha, j = 1)
    alpha_star_2 <- alpha_star_j(rho = rho, alpha = alpha, j = 2)

    tau_star_1 <- tau_star_j(rho = rho, alpha = alpha, nu = nu, j = 1)
    tau_star_2 <- tau_star_j(rho = rho, alpha = alpha, nu = nu, j = 2)

    I1 <- pest(x = comp1_1, scale = 1, shape = alpha_star_1, extended = tau_star_1, df = nu + 1)
    I2 <- pest(x = comp1_2, scale = 1, shape = alpha_star_2, extended = tau_star_2, df = nu + 1)

    if (type == "lower") {
      return(exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    } else if (type == "inv.lower") {
      return(1 - exp(-1 / q[1] * I1 - 1 / q[2] * I2))
    } else if (type == "upper") {
      return(1 / q[1] * (1 - I1) + 1 / q[2] * (1 - I2))
    }
  }

  p.est.3d <- function(q, par, type) {
    dim <- length(q)
    if (dim != 3) {
      stop("Dimension should be 3")
    }

    rho <- par[1:3]
    alpha <- par[4:6]
    nu <- par[7]

    if ((any(rho <= -1) || any(rho >= 1) || (nu <= 0)) == 1) {
      stop("Extremal-t parameters ill defined")
    }

    Sigma_bar_1 <- Sigma_bar_j(rho = rho, j = 1)
    Sigma_bar_2 <- Sigma_bar_j(rho = rho, j = 2)
    Sigma_bar_3 <- Sigma_bar_j(rho = rho, j = 3)

    if (sum(eigen(Sigma_bar_1)$values < 0) >= 1) {
      stop("negative definite partial correlation matrix")
    }
    if (sum(eigen(Sigma_bar_2)$values < 0) >= 1) {
      stop("negative definite partial correlation matrix")
    }
    if (sum(eigen(Sigma_bar_3)$values < 0) >= 1) {
      stop("negative definite partial correlation matrix")
    }

    nu_1 <- nu_p(rho = rho, alpha = alpha, nu = nu, j = 1)
    nu_2 <- nu_p(rho = rho, alpha = alpha, nu = nu, j = 2)
    nu_3 <- nu_p(rho = rho, alpha = alpha, nu = nu, j = 3)

    x1_bar <- q[1] * nu_1
    x2_bar <- q[2] * nu_2
    x3_bar <- q[3] * nu_3

    comp1_1 <- sqrt((nu + 1) / (1 - rho[1]^2)) * ((x2_bar / x1_bar)^(1 / nu) - rho[1])
    comp2_1 <- sqrt((nu + 1) / (1 - rho[2]^2)) * ((x3_bar / x1_bar)^(1 / nu) - rho[2])
    comp1_2 <- sqrt((nu + 1) / (1 - rho[1]^2)) * ((x1_bar / x2_bar)^(1 / nu) - rho[1])
    comp2_2 <- sqrt((nu + 1) / (1 - rho[3]^2)) * ((x3_bar / x2_bar)^(1 / nu) - rho[3])
    comp1_3 <- sqrt((nu + 1) / (1 - rho[2]^2)) * ((x1_bar / x3_bar)^(1 / nu) - rho[2])
    comp2_3 <- sqrt((nu + 1) / (1 - rho[3]^2)) * ((x2_bar / x3_bar)^(1 / nu) - rho[3])

    alpha_star_1 <- alpha_star_j(rho = rho, alpha = alpha, j = 1)
    alpha_star_2 <- alpha_star_j(rho = rho, alpha = alpha, j = 2)
    alpha_star_3 <- alpha_star_j(rho = rho, alpha = alpha, j = 3)

    tau_star_1 <- tau_star_j(rho = rho, alpha = alpha, nu = nu, j = 1)
    tau_star_2 <- tau_star_j(rho = rho, alpha = alpha, nu = nu, j = 2)
    tau_star_3 <- tau_star_j(rho = rho, alpha = alpha, nu = nu, j = 3)

    I1 <- pmest(x = c(comp1_1, comp2_1), scale = Sigma_bar_1, shape = alpha_star_1, extended = tau_star_1, df = nu + 1)
    I2 <- pmest(x = c(comp1_2, comp2_2), scale = Sigma_bar_2, shape = alpha_star_2, extended = tau_star_2, df = nu + 1)
    I3 <- pmest(x = c(comp1_3, comp2_3), scale = Sigma_bar_3, shape = alpha_star_3, extended = tau_star_3, df = nu + 1)

    if (type == "lower") {
      return(exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    } else if (type == "inv.lower") {
      return(1 - exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    } else if (type == "upper") {
      return(1 - exp(-1 / q[1]) - exp(-1 / q[2]) - exp(-1 / q[3])
        + p.est.2d(q[1:2], c(rho[1], alpha[1:2], nu), type = "lower") + p.est.2d(q[c(1, 3)], c(rho[2], alpha[c(1, 3)], nu), type = "lower")
        + p.est.2d(q[2:3], c(rho[3], alpha[2:3], nu), type = "lower") - exp(-1 / q[1] * I1 - 1 / q[2] * I2 - 1 / q[3] * I3))
    }
  }

  dim <- length(q)
  if (!dim_ExtDep(model = "EST", par = par, dim = dim)) {
    stop("Length of par is incorrect")
  }
  if (dim == 2) {
    return(p.est.2d(q = q, par = par, type = type))
  } else if (dim == 3) {
    return(p.est.3d(q = q, par = par, type = type))
  }
}

### Definition of the angular probability measure

ph <- function(w, beta) {
  k <- length(beta) - 1
  j <- 1:k
  #  res <- diff(beta) * dbeta(w, j, k-j+1)
  res <- 0.5 * (diff(beta) + 1 / k) * dbeta(w, j, k - j + 1)
  return(sum(res))
}


#########################################
#########################################
### END Internal functions for pExtDep
#########################################
#########################################