# Joe's Nu-Skew and the copBasic Nu-Star of a Copula 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 [REMOVE_ME]$$ \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{.} [REMOVE_ME_2]$$ 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 [REMOVE_ME]$$ \nu^\star_\mathbf{C} = 12\int\!\!\int_{\mathcal{I}^2} (v+u)\mathbf{C}(u,v)\, \mathrm{d}u\mathrm{d}v - 4\mbox{,} [REMOVE_ME_2]$$ 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). utf8 ## 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). ```r 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 `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. ## Returns The value for $\nu_\mathbf{C}$ or $\nu^\star_\mathbf{C}$ is returned. ## 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 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) ``` ## 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 . ## References Joe, H., 2014, Dependence modeling with copulas: Boca Raton, CRC Press, 462 p. ## Author(s) W.H. Asquith ## See Also `uvskew`, `blomCOP`, `footCOP`, `giniCOP`, `hoefCOP`, `rhoCOP`, `tauCOP`, `wolfCOP` ## Examples ```r 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) ```

joeskewCOP function