joeskewCOP: Joe's Nu-Skew and the copBasic Nu-Star of a Copula

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

Joe's Nu-Skew and the copBasic Nu-Star of a Copula

Description

Compute the measure of permutation asymmetry, which can be thought of as bivariate skewness, named for the copBasic package as Nu-Skew \nu_\mathbf{C} of a copula according to Joe (2014, p. 66) by

\nu_\mathbf{C} = 3\mathrm{E}[UV^2 - U^2V] = 6\int\!\!\int_{\mathcal{I}^2} (v-u)\mathbf{C}(u,v)\, \mathrm{d}u\mathrm{d}v\mbox{.}

This definition is effectively the type="nu" for the function for which the multiplier 6 has been converted to 96 as explained in the Note.

Numerical results indicate \nu_\mathbf{W} \approx 0 (W), \nu_\mathbf{\Pi} = 0 (P), \nu_\mathbf{M} \approx 0 (M), \nu_\mathbf{PL} \approx 0 for all \Theta (PLcop), and the \nu^\star_\mathbf{GH} = 0 (GHcop); copulas with mirror symmetry across the equal value line have \nu_\mathbf{C} = 0.

Asymmetric copulas do exist. For example, consider an asymmetric Gumbel–Hougaard \mathbf{GH} copula with \Theta_p = (5,0.8,p):

  optimize(function(p) { nuskewCOP(cop=GHcop, para=c(5,0.8, p)) },
           c(0,0.99) )$minimum
  UV <- simCOP(n=10000, cop=GHcop, c(5,0.8, 0.2836485)) # inspect the graphics
  48*mean(UV$U*$V^2 - UV$U^2*UV$V) # -0.2847953 (not the 3rd parameter)

The minimization yields \nu_{\mathbf{GH}(5, 0.8, 0.2836485)} = -0.2796104, which is close the expectation computed where 48 = 96/2.

A complementary definition is supported, triggered by type="nustar", and is computed by

\nu^\star_\mathbf{C} = 12\int\!\!\int_{\mathcal{I}^2} (v+u)\mathbf{C}(u,v)\, \mathrm{d}u\mathrm{d}v - 4\mbox{,}

which has been for the copBasic package, \nu^\star_\mathbf{C} is named as Nu-Star, which the 12 and the -4 have been chosen so that numerical results indicate \nu^\star_\mathbf{W} = -1 (W), \nu^\star_\mathbf{\Pi} = 0 (P), and \nu^\star_\mathbf{M} = +1 (M).

Lastly, the uvlmoms function provides for a quantile-based measure of bivariate skewness based on the difference U - V that also is discussed by Joe (2014, p. 66).

Usage

joeskewCOP(cop=NULL, para=NULL, type=c("nu", "nustar", "nuskew"),
                               as.sample=FALSE, brute=FALSE, delta=0.002, ...)

nuskewCOP(cop=NULL, para=NULL, as.sample=FALSE, brute=FALSE, delta=0.002, ...)
nustarCOP(cop=NULL, para=NULL, as.sample=FALSE, brute=FALSE, delta=0.002, ...)

Arguments

cop	A copula function;
para	Vector of parameters or other data structure, if needed, to pass to the copula;
type	The type of metric to compute (nu and nuskew are synonymous for \nu_\mathbf{C} and nustar is for \nu^\star_\mathbf{C});
brute	Should brute force be used instead of two nested integrate() functions to perform the double integration;
delta	The \mathrm{d}u and \mathrm{d}v for the brute force integration using brute;
as.sample	A logical controlling whether an optional R data.frame in para is used to compute the sample \hat\nu or \hat\nu^\star (see Note). If set to -1, then the message concerning CPU effort will be surpressed; and
...	Additional arguments to pass.

Details

The implementation of joeskewCOP for copBasic provides the second metric of asymmetry, but why? Consider the results that follow:

  joeskewCOP(cop=GHcop, para=c(5, 0.8,    0.2836485), type="nu")
     # -0.2796104
  joeskewCOP(cop=GHcop, para=c(5, 0.2836485,    0.8), type="nu")
     # +0.2796103
  joeskewCOP(cop=GHcop, para=c(5, 0.8,    0.2836485), type="nu")
     #  0.3571276
  joeskewCOP(cop=GHcop, para=c(5, 0.2836485,    0.8), type="nu")
     #  0.3571279
  tauCOP(    cop=GHcop, para=c(5, 0.2836485,    0.8))
     #  0.2443377

The demonstration shows—at least for the symmetry (switchability) of the 2nd and 3rd parameters (\pi_2 and \pi_3) of the asymmetric \mathbf{GH} copula—that the first definition \nu is magnitude symmetric but carries a sign change. The demonstration shows magnitude and sign stability for \nu^\star, and ends with Kendall Tau (tauCOP). Collectively, Kendall Tau (or the other symmetric measures of association, e.g. blomCOP, footCOP, giniCOP, hoefCOP, rhoCOP, wolfCOP) when combined with \nu and \nu^\star might provide a framework for parameter optimization of the asymmetric \mathbf{GH} copula (see below).

The asymmetric \mathbf{GH}_{(5, 0.2836485, 0.8)} is not radial (isCOP.radsym) or permutation (isCOP.permsym), but if \pi_2 = \pi_3 then the resulting \mathbf{GH} copula is not radially symmetric but is permutation symmetric:

  isCOP.radsym( cop=GHcop, para=c(5, 0.2836485, 0.8)) # FALSE
  isCOP.permsym(cop=GHcop, para=c(5, 0.2836485, 0.8)) # FALSE
  isCOP.radsym( cop=GHcop, para=c(5, 0.8,       0.8)) # FALSE
  isCOP.permsym(cop=GHcop, para=c(5, 0.8,       0.8)) # TRUE

The use of \nu_\mathbf{C} and \nu^\star_\mathbf{C} with a measure of association is just suggested above for parameter optimization. Suppose we have \mathbf{GH}_{(5,0.5,0.7)} with Spearman Rho \rho = 0.4888, \nu = 0.001475, and \nu^\star = 0.04223, and the asymmetric \mathbf{GH} coupla is to be fit. Parameter estimation for the asymmetric \mathbf{GH} is accomplished by

  "fitGHcop" <- function(hats, assocfunc=rhoCOP, init=NA, eps=1E-4, ...) {
     H <- GHcop # shorthand for the copula
     "objfunc" <- function(par) {
        par[1]   <- ifelse(par[1] < 1, return(Inf), exp(par[1])) # edge check
        par[2:3] <-  pnorm(par[2:3]) # detransform
        hp <- c(assocfunc(H, par), nuskewCOP(H, par), nustarCOP(H, par))
        return(sum((hats-hp)^2))
     }
     # Theta=1 and Pi2 = Pi3 = 1/2 # as default initial estimates
     if(is.na(init)) init <- c(1, rep(1/2, times=2))
     opt  <- optim(init, objfunc, ...); par <- opt$par
     para <- c( exp(par[1]), pnorm(par[2:3]) )
     names(para) <- c("Theta", "Pi2", "Pi3")
     fit <- c(assocfunc(H, para), nuskewCOP(H, para), nustarCOP(H, para))
     txt <- c("AssocMeasure", "NuSkew", "NuStar")
     names(fit) <- txt; names(hats) <- txt
     if(opt$value > eps) warning("inspect the fit")
     return(list(para=para, fit=fit, given=hats, optim=opt))
  }
  father <- c(5,.5,.7)
  densityCOPplot(cop=GHcop, para=father, contour.col=8)
  fRho  <- rhoCOP(   cop=GHcop, father)
  fNu   <- nuskewCOP(cop=GHcop, father)
  fStar <- nustarCOP(cop=GHcop, father)

  child <- fitGHcop(c(fRho, fNu, fStar))$para
  densityCOPplot(cop=GHcop, para=child, ploton=FALSE)

  cRho  <- rhoCOP(   cop=GHcop, child)
  cNu   <- nuskewCOP(cop=GHcop, child)
  cStar <- nustarCOP(cop=GHcop, child)
  message("Father stats: ", paste(fRho, fNu, fStar, sep=", "))
  message("Child  stats: ", paste(cRho, cNu, cStar, sep=", "))
  message("Father para: ",  paste(father,      collapse=", "))
  message("Child  para: ",  paste(child,       collapse=", "))

The initial parameter estimate has the value \Theta = 1, which is independence for the one parameter \mathbf{GH}. The two other parameters are set as \pi_2 = \pi_3 = 1/2 to be in the mid-point of their domain. The transformations using the log() \leftrightarrow exp() and qnorm() \leftrightarrow pnorm() functions in R are used to keep the optimization in the viable parameter domain. The results produce a fitted copula of \mathbf{GH}_{(4.907, 0.5006, 0.7014)}. This fit aligns well with the parent, and the three statistics are essentially matched during the fitting.

The \nu^\star_\mathbf{C} can be similar to rhoCOP, but differences do exist. In the presence of radial symmetry, (\nu_\mathbf{C} == 0), the \nu^\star_\mathbf{C} is nearly equal to Spearman Rho for some copulas. Let us test further:

  p <- 10^seq(0,2,by=.01)
  s <- sapply(p, function(t) nustarCOP(cop=GHcop, para=c(t)))
  r <- sapply(p, function(t)    rhoCOP(cop=GHcop, para=c(t)))
  plot(p,s, log="x", type="l", col=3, lwd=3); lines(p,r)

Now let us add some asymmetry

  s <- sapply(p, function(t) nustarCOP(cop=GHcop, para=c(t, 0.25, 0.75)))
  r <- sapply(p, function(t)    rhoCOP(cop=GHcop, para=c(t, 0.25, 0.75)))
  plot(p,s, log="x", type="l", col=3, lwd=3); lines(p,r)

Now let us choose a different (the Clayton) copula

  s <- sapply(p, function(t) nustarCOP(cop=CLcop, para=c(t)))
  r <- sapply(p, function(t)    rhoCOP(cop=CLcop, para=c(t)))
  plot(p,s, log="x", type="l", col=3, lwd=3); lines(p,r)

Value

The value for \nu_\mathbf{C} or \nu^\star_\mathbf{C} is returned.

Note

The \nu_\mathbf{C} definition is given with a multiplier of 6 on the integrals in order to agree with Joe (2014) relation that is also shown. However, in mutual parameter estimation experiments using a simple sum-of-square errors as shown in the Details, it is preferred to have \nu_\mathbf{C} measured on a larger scale. Where does the 96 then come from? It is heuristically made so that the upright and rotated cophalf (see Examples under asCOP and bilmoms for this copula) acquires \nu_\mathbf{C} values of +1 and -1, respectively. As a result to make back comparisons to Joe results, the ratios of 96 are made in this documentation.

The source code shows slightly different styles of division by the sample size as part of the sample estimation of the statistics. The \hat\nu using just division by the sample size as testing indicates that this statistic is reasonably unbiased for simple copula. The \hat\nu^\star with similar division is a biased statistic and the bias is not symmetrical in magnitude or sign it seems whether the \hat\nu^\star is positive or negative. The salient code is spm <- ifelse(corsgn == -1, +2.4, +1.1) within the sources for which the corrections were determined heuristically through simulation, and corsgn is the sign of the sample Spearman Rho through the cor() function of R.

Author(s)

W.H. Asquith

References

Joe, H., 2014, Dependence modeling with copulas: Boca Raton, CRC Press, 462 p.

See Also

uvskew, blomCOP, footCOP, giniCOP, hoefCOP, rhoCOP, tauCOP, wolfCOP

Examples

nuskewCOP(cop=GHcop,para=c(1.43,1/2,1))*(6/96) # 0.005886 (Joe, 2014, p. 184; 0.0059)

## Not run: 
joeskewCOP(            cop=GHcop, para=c(8, .7, .5)) # -0.1523491
joeskewCOP(            cop=GHcop, para=c(8, .5, .7)) # +0.1523472
# UV <- simCOP(n=1000, cop=GHcop, para=c(8, .7, .5)) # see the switch in
# UV <- simCOP(n=1000, cop=GHcop, para=c(8, .5, .7)) # curvature
## End(Not run)

## Not run: 
para=c(19,0.3,0.8); set.seed(341)
nuskew <-  nuskewCOP( cop=GHcop, para=para) # 0.3057744
UV <- simCOP(n=10000, cop=GHcop, para=para) #   a large simulation
mean((UV$U - UV$V)^3)/(6/96)                # 0.3127398

# Two other definitions of skewness follow and are not numerically the same.
uvskew(u=UV$U, v=UV$V, umv=TRUE)  # 0.3738987  (see documentation uvskew)
uvskew(u=UV$U, v=UV$V, umv=FALSE) # 0.3592739  ( or documentation uvlmoms)
# Yet another definition of skew, which requires large sample approximation
# using the L-comoments (3rd L-comoment is L-coskew).
lmomco::lcomoms2(UV)$T3 # L-coskew of the simulated values [1,2] and [2,1]
#             [,1]        [,2]
#[1,]  0.007398438  0.17076600
#[2,] -0.061060260 -0.00006613
# See the asymmetry in the two L-coskew values and consider this in light of
# the graphic produced by the simCOP() called for n=10,000. The T3[1,1] is
# the sampled L-skew (univariate) of the U margin and T3[2,2] is the same
# but for the V margin. Because the margins are uniform (ideally) then these
# for suitable large sample must be zero because the L-skew of the uniform
# distribution is by definition zero.
#
# Now let us check the sample estimator for sample of size n=300, and the
# t-test will (should) result in acceptance of the NULL hypothesis.
S <- replicate(60, nuskewCOP(para=simCOP(n=300, cop=GHcop, para=para,
                                         graphics=FALSE), as.sample=TRUE))
t.test(S, mu=nuskew)
# t = 0.004633, df = 59, p-value = 0.9963
# alternative hypothesis: true mean is not equal to 0.3057744
# 95 percent confidence interval:
#  0.2854282 0.3262150
# sample estimates:
# mean of x
# 0.3058216 
## End(Not run)

## Not run: 
# Let us run a large ensemble of copula properties that use the whole copula
# (not tail properties). We composite a Plackett with a Gumbel-Hougaard for
# which the over all association (correlation) sign is negative, but amongst
# these statistics with nuskew and nustar at the bottom, there are various
# quantities that can be extracted. These could be used for fitting.
set.seed(873)
para <- list(cop1=PLcop, cop2=GHcop, alpha=0.6, beta=0.9,
             para1=.005, para2=c(8.3,0.25,0.7))
UV <- simCOP(1000, cop=composite2COP, para=para) # just to show
  blomCOP(composite2COP, para)            # -0.4078657
  footCOP(composite2COP, para)            # -0.2854227
  hoefCOP(composite2COP, para)            # +0.5713775
  lcomCOP(composite2COP, para)$lcomUV[3]  # +0.1816084
  lcomCOP(composite2COP, para)$lcomVU[3]  # +0.1279844
   rhoCOP(composite2COP, para)            # -0.5688417
rhobevCOP(composite2COP, para)            # -0.2005210
   tauCOP(composite2COP, para)            # -0.4514693
  wolfCOP(composite2COP, para)            # +0.5691933
nustarCOP(composite2COP, para)            # -0.5172434
nuskewCOP(composite2COP, para)            # +0.0714987 
## End(Not run)

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