joint.curvesCOP: Compute Coordinates of the Marginal Probabilities given joint...

In wasquith/copBasic: General Bivariate Copula Theory and Many Utility Functions

Compute Coordinates of the Marginal Probabilities given joint AND or OR Probabilities

Description

Compute the coordinates of the bivariate marginal probabilities for variables U and V given selected probabilities levels t for a copula \mathbf{C}(u,v) for v with respect to u. For the case of a joint and probability, symbolically the solution is

\mathrm{Pr}[U \le v,\ V \le v] = t = \mathbf{C}(u,v)\mbox{,}

where U \mapsto [t_i, u_{j}, u_{j+1}, \cdots, 1; \Delta t] (an irregular sequence of u values from the ith value of t_i provided through to unity) and thus

t_i \mapsto \mathbf{C}(u=U, v)\mbox{,}

and solving for the sequence of v. The index j is to indicate that a separate loop is involved and is distinct from i. The pairings \{u(t_i), v(t_i)\} for each t are packaged as an R data.frame. This operation is very similiar to the plotting capabilities in level.curvesCOP for level curves (Nelsen, 2006, pp. 12–13) but implemented in the function joint.curvesCOP for alternative utility.

For the case of a joint or probability, the dual of a copula (function) or \tilde{\mathbf{C}}(u,v) from a copula (Nelsen, 2006, pp. 33–34) is used and symbolically the solution is:

\mathrm{Pr}[U \le v \mathrm{\ or\ } V \le v] = t = \tilde{\mathbf{C}}(u,v) = u + v - \mathbf{C}(u,v)\mbox{,}

where U \mapsto [0, u_j, u_{j+1}, \cdots, t_i; \Delta t] (an irregular sequence of u values from zero through to the ith value of t) and thus

t_i \mapsto \tilde{\mathbf{C}}(u=U, v)\mbox{,}

and solving for the sequence of v. The index j is to indicate that a separate loop is involved and is distinct from i. The pairings \{u(t_i), v(t_i)\} for each t are packaged as an R data.frame.

Usage

joint.curvesCOP(cop=NULL, para=NULL, type=c("and", "or"),
                probs=c(0.5, 0.8, 0.90, 0.96, 0.98, 0.99, 0.995, 0.998),
                zero2small=TRUE, small=1E-6, divisor=100, delu=0.001, ...)

Arguments

cop	A copula function;
para	Vector of parameters or other data structure, if needed, to pass to the copula;
type	What type of joint probability is to be computed;
probs	The joint probabilities t_i from which to compute the coordinates. The default values represent especially useful annual return period equivalents that are useful in hydrologic risk analyses;
zero2small	A logical controlling whether exactly zero value for probability are converted to a small value and exactly unity values for probability are converted to the value 1 - small; this logical is useful if transformation from probability space into standard normal variates or Gumbel reduced variates (see function prob2grv() in package lmomco) is later desired by the user for attendant graphics (see Examples section);
small	The value for small described for zero2small;
divisor	A divisor on a computation of a \Delta t for incrementing through the irregularly-spaced u domain as part of the coordinate computation (see source code);
delu	A \Delta u for setup of the incrementing through the irregularly-space u domain as part of the coordinate computation (see source code); and
...	Additional arguments to pass to the duCOP function of copBasic or uniroot() function in R.

Value

An R list is returned with elements each of the given probs.

Note

The arguments divisor and delu provide flexibility to obtain sufficient smoothness in the coordinate curvatures for a given t. The pairings \{u(t_i), v(t_i)\} for each t packaged as data.fames within the returned list each have their own unique length, and this is the reason that a single “master” data.frame is not returned by this function.

Extended Example—The code below shows both types of joint probability being computed using the default probs. The plotting is made in Gumbel reduced variates (GRV; see prob2grv in package lmomco). This transformation is somewhat suitable for the magnitude variation in and at tail depth of the probs. Also with transformation is being used, the zero2small logical is kept at TRUE, which is appropriate. The zero2small being set is also useful if standard normal variate transformation (by the qnorm() function in R) were used instead.

The Gumbel–Hougaard copula (GHcop) and a reversed Gumbel–Hougaard copula rGH are composited together by composite2COP. These were chosen so that some asymmetry in the solutions by joint.curvesCOP could be seen. We begin by specifying symmetrical plotting limits for later use and then creating a function for the reversed Gumbel–Hougaard and setting the parameters and composition weights:

  grvlim <- lmomco::prob2grv(c(0.25,0.999)) # out to 1,000 years
  "rGHcop" <- function(u,v, ...) { u + v - 1 + GHcop(1-u, 1-v, ...) }
  para <- list(alpha=0.16, beta=0.67, cop1 =GHcop, cop2 =rGHcop,
                                      para1=4.5,   para2=2.2)
  tau    <- tauCOP(    cop=composite2COP, para=para) # Tau    = 0.351219
  nuskew <- nuskewCOP(cop=composite2COP, para=para)  # Nuskew = 0.084262
  UV <- simCOP(n=1000,  cop=composite2COP, para=para, snv=TRUE)

The code also computed the Kendall Tau (tauCOP) and Nu-Skew (nuskewCOP) and the results shown. The code finishes with a simulation by simCOP of the copula composition just for reference.

Next, we compute and plot the joint probability curves. The tolerance for the uniroot calls is reset from the R defaults because slight wooble in the numerical solutions exists otherwise. The AND and OR lists each provide data.frames from which successive graphic calls plot the lines. The second plot() call is commented out so that both sets of joint probability curves are drawn on the same plot.

  AND <- joint.curvesCOP(cop=composite2COP, para=para, type="and",
                         divisor=1000, tol=.Machine$double.eps)
  plot(grvlim, grvlim, type="n",
       xlab="GUMBEL REDUCED VARIATE IN U", ylab="GUMBEL REDUCED VARIATE IN V")
  for(t in sort(as.numeric(names(AND)))) {
      UV <- get(as.character(t), AND)
      lines(lmomco::prob2grv(UV$U),         lmomco::prob2grv(       UV$V))
      text( lmomco::prob2grv(median(UV$U)), lmomco::prob2grv(median(UV$V)),
           as.character(round(1/(1-t)), digits=0))
  }

  OR <- joint.curvesCOP(cop=composite2COP, para=para, type="or",
                        divisor=1000, tol=.Machine$double.eps)
  for(t in sort(as.numeric(names(OR)))) {
     UV <- get(as.character(t), OR)
     lines(lmomco::prob2grv(UV$U), lmomco::prob2grv(UV$V), col=2)
     text( lmomco::prob2grv(median(UV$U)), lmomco::prob2grv(median(UV$V)),
          as.character(round(1/(1-t)), digits=0), col=2)
  }
  mtext("Return Periods: black=cooperative risk, red=dual risk")
  abline(0,1, lty=2) # dash line is simply and equal value line

Black Curves—The black curves represent the nonexceedance joint and condition. The black curves are a form of level curve (see level.curvesCOP), but it seems appropriate to not name them as such because Nelsen's examples and others usually have the level curves on an even step interval of probability such as 10-percent level curves. The complement of nonexceedance joint and is the probability level that either or both random variables (say “hazards”) U or V causes “failure” at the respective return period.

Red Curves—The red curves represent the nonexceedance joint or (inclusive) condition. The complement of nonexceedance joint or (inclusive) is the probability level that both random variables (say “hazards”) U or V must simultaneously (or dually) occur for “failure” at the respective return period. Note, there is obviously asymmetry in the joint or curves.

Interpretation—Because the two hazards can “cooperate” to cause failure (see coCOP) for an equal level of protection (say 500-year event) relative to the complement of nonexceedance joint or (inclusive) condition (see surCOP), the marginal probabilities must be considerably higher. Users are encouraged to review copBasic-package and the figure therein.

Author(s)

W.H. Asquith

References

Nelsen, R.B., 2006, An introduction to copulas: New York, Springer, 269 p.

See Also

diagCOPatf, duCOP, jointCOP, joint.curvesCOP2, level.curvesCOP

Examples

# See Note section

wasquith/copBasic documentation built on Dec. 13, 2024, 6:39 p.m.