bilmoms {copBasic}

Bivariate L-moments and L-comoments of a Copula

Description

Attention: This function is deprecated in favor of lcomCOP, which uses only direct numerical integrate() on the integrals shown below. The bilmoms function is strictly based on Monte Carlo integration.

Compute the bivariate L-moments (ratios) (\delta^{[\ldots]}_{k;\mathbf{C}}) of a copula \mathbf{C}(u,v; \Theta) and remap these into the L-comoment matrix counterparts (Serfling and Xiao, 2007; Asquith, 2011) including L-correlation (Spearman Rho), L-coskew, and L-cokurtosis. As described by Brahimi et al. (2015), the first four bivariate L-moments \delta^{[12]}_k for random variable X^{(1)} or U with respect to (wrt) random variable X^{(2)} or V are defined as

\delta^{[12]}_{1;\mathbf{C}} = 2\int\!\!\int_{\mathcal{I}^2} \mathbf{C}(u,v)\,\mathrm{d}u\mathrm{d}v - \frac{1}{2}\mbox{,}

\delta^{[12]}_{2;\mathbf{C}} = \int\!\!\int_{\mathcal{I}^2} (12v - 6) \mathbf{C}(u,v)\,\mathrm{d}u\mathrm{d}v - \frac{1}{2}\mbox{,}

\delta^{[12]}_{3;\mathbf{C}} = \int\!\!\int_{\mathcal{I}^2} (60v^2 - 60v + 12) \mathbf{C}(u,v)\,\mathrm{d}u\mathrm{d}v - \frac{1}{2}\mbox{, and}

\delta^{[12]}_{4;\mathbf{C}} = \int\!\!\int_{\mathcal{I}^2} (280v^3 - 420v^2 + 180v - 20) \mathbf{C}(u,v)\,\mathrm{d}u\mathrm{d}v - \frac{1}{2}\mbox{,}

where the bivariate L-moments are related to the L-comoment ratios by

6\delta^{[12]}_k = \tau^{[12]}_{k+1}\mbox{\quad and \quad}6\delta^{[21]}_k = \tau^{[21]}_{k+1}\mbox{,}

where in otherwords, “the third bivariate L-moment \delta^{[12]}_3 is one sixth the L-cokurtosis \tau^{[12]}_4.” The first four bivariate L-moments yield the first five L-comoments (there is no first order L-comoment ratio). The terms and nomenclature are not easy and also the English grammar adjective “ratios” is not always consistent in the literature. The \delta^{[\ldots]}_{k;\mathbf{C}} are ratios, and the returned bilcomoms element by this function holds matrices for the marginal means, marginal L-scales and L-coscales, and then the ratio L-comoments.

Similarly, the \delta^{[21]}_k are computed by switching u \rightarrow v in the polynomials within the above integrals multiplied to the copula in the system of equations with u. In general, \delta^{[12]}_k \not= \delta^{[21]}_k for k > 1 unless in the case of permutation symmetric (isCOP.permsym) copulas. By theory, \delta^{[12]}_1 = \delta^{[21]}_1 = \rho_\mathbf{C}/6 where \rho_\mathbf{C} is a Spearman Rho rhoCOP.

The integral for \delta^{[12]}_{4;\mathbf{C}} does not appear in Brahimi et al. (2015) but this and the other forms are verified in the Examples and discussion in Note. The four k \in [1,2,3,4] for U wrt V and V wrt U comprise a full spectrum of system of seven (not eight) equations. One equation is lost because \delta^{[12]}_1 = \delta^{[21]}_1.

Usage

bilmoms(cop=NULL, para=NULL, only.bilmoms=FALSE, n=1E5,
                  sobol=TRUE, scrambling=0, ...)

Arguments

cop

A copula function;

para

Vector of parameters or other data structure, if needed, to pass to the copula;

only.bilmoms

A logical to trigger return of the \delta_k and skip L-comoment computation;

n

The Monte Carlo integration size. The default seems to be at least an order of magnitude greater than needed for many applied problems;

sobol

A logical trigging Sobol sequences for the Monte Carlo integration instead of the bivariate uniform distribution (independence). The Sobol sequences are dependent on the sobol() function of the randtoolbox package, and the Sobol sequences canvas the \mathcal{I}^2 domain for smaller n values than required if statistical independence is used for the Monte Carlo integration. Note, randtoolbox at least at version 2.0.+ has “scrambling” of Sobol sequences temporarily disabled, and hence scrambling=0 as default for bilmoms;

scrambling

The argument of the same name for randtoolbox::sobol(); and

...

Additional arguments to pass to the densityCOP function.

Value

An R list of the bivariate L-moments is returned.

bilmomUV

The bivariate L-moments \delta^{[12]}_k of U with respect to V for k \in [1,2,3,4];

bilmomVU

The bivariate L-moments \delta^{[21]}_k of V with respect to U for k \in [1,2,3,4];

error.rho

An “error” term in units of \delta^{[12 \& 21]}_1 used to judge whether the sample size for the Monte Carlo integration is sufficient based on a comparison to the Spearman Rho from direct numerical integration (not Monte Carlo based) using rhoCOP of the copula. Values for error.rho < 1E{-}5 seem to be sufficient to judge whether n is large enough;

bilcomoms

If not only.bilmoms, another R list holding the L-comoments (see Note) computed by simple remapping of the \delta^{[\ldots]}_k and parallel in structure to the function lcomoms2() of the lmomco package; and

source

An attribute identifying the computational source of the bivariate L-moments and bivariate L-comoments: “bilmoms.”

Note

The mapping of the bivariate L-moments to their L-comoment matrix counterparts is simple but nuances should be discussed and the meaning of the error.rho needs further description. The extra effort to form L-comoment matrices (Serfling and Xiao, 2007; Asquith, 2011) is made so that output matches the structure of the sample L-comoment matrices from the lcomoms2() function of the lmomco package.

Concerning the triangular or tent-shaped copula of Nelsen (2006, exer. 3.7, pp. 64–65) for demonstration, simulate from the triangular copula a sample of size m = 20{,}000 and compute some sample L-comoments using the following CPU intensive code. The function asCOP completes the vectorization needed for non-Monte Carlo integration for rhoCOP.

  "trianglecop" <- function(u,v, para=NULL, ...) {
    # If para is set, then the triangle is rotated 90d clockwise.
    if(! is.null(para) && para == 1) { t <- u; u <- v; v <- t }
    if(length(u) > 1 | length(v) > 1) stop("only scalars for this function")
    v2<-v/2; if(0   <= u    & u    <= v2 & v2 <= 1/2) { return(u    )
    } else   if(0   <=   v2 &   v2 <  u  &  u < 1-v2) { return(v2   )
    } else   if(1/2 <= 1-v2 & 1-v2 <= u  &  u <= 1  ) { return(u+v-1)
    } else { stop("should not be here in logic") }
  }
  "TriCop" <- function(u,v, ...) { asCOP(u,v, f=trianglecop, ...) }
  m=20000; SampleUV <- simCOP(n=m, cop=TriCop, graphics=FALSE)
  samLC  <- lcomoms2(SampleUV, nmom=5)
  theoLC <- bilmoms(cop=TriCop)

The Error in Rho Computation—The \rho_\mathbf{C} of the copula by numerical integration is computed internally to bilmoms as

  rhoC <- rhoCOP(cop=TriCop) # -1.733858e-17

and used to compute the error.rho for bilmoms (see next code snippet). The \rho_\mathbf{C} is obviously zero for this copula. Therefore, the bivariate association of TriCop is zero and thus is an example of a perfectly dependent situation yet of zero correlation. The bivariate L-moments and L-comoments of this copula are computed as

  mean(replicate(20, bilmoms(cop=TriCop)$error.rho)) # 7.650723e-06

where the error.rho is repeated trials appears firmly <1e-5, which is near zero (\epsilon_\rho \approx 0). The error.rho term is defined by taking the first bivariate L-moment and numerically integrated \rho_\mathbf{C} through rhoCOP and computing the terms

\epsilon^{[12]}_\rho = |\delta^{[12]}_{1;\mathbf{C}} - (\rho_\mathbf{C}/6)|\mbox{,}

\epsilon^{[21]}_\rho = |\delta^{[21]}_{1;\mathbf{C}} - (\rho_\mathbf{C}/6)|\mbox{, and}

\epsilon_\rho = \frac{\epsilon^{[12]}_\rho + \epsilon^{[21]}_\rho}{2}\mbox{,}

where the error.rho = \epsilon_\rho, and values near zero are obviously favorable because this indicates that the Monte Carlo integration sample size n argument is sufficiently large to effectively canvas the \mathcal{I}^2 domain. For the situation here, the theoretical \rho_\mathbf{C} = 0, but for n = n = 100, the error.rho \approx 0.006 (e.g. bilmoms(n=100, cop=TriCop)$error.rho) through a 20-unit replication, which is a hint that 100 samples are not large enough and that should be obvious.

The reasoning behind using the error.rho between conventional numerical integration and the Monte Carlo integration (error.rho) is that \rho_\mathbf{C} is symmetrical. This choice of “convergence” assessment reduces somewhat the sample size needed for Monte Carlo integration into single number representing error.

Discussion of Theoretical L-comoments—The theoretical L-comoments in the format structure of the sample L-comoments by the lcomoms2() function of the lmomco package are formed by the bilmoms function, and the theoretical values are shown below in sequence with details listed by L-comoment. Now we extract the L-comoment matrices and show the first L-comoment matrix (the matrix of the means):

  theoLClcm <- theoLC$bilcomoms
  print(theoLClcm$L1)
            [,1]      [,2]
  [1,] 0.4999939        NA
  [2,]        NA 0.5000032

where the diagonal should be filled with 1/2, if the n is suitably large, because 1/2 is the mean of the marginal uniform random variables. By definition the secondary diagonal has NAs. The values shown above are extremely close supporting the idea that default n is large enough. The matrix of means is otherwise uninformative.

The second L-comoment matrix (L-scales and L-coscales) is

  print(theoLClcm$L2)
                [,1]          [,2]
  [1,]  1.666677e-01 -3.466202e-06
  [2,] -3.466209e-06  1.666674e-01

where the diagonal should be filled with 1/6, if the n is suitably large, because 1/6 is the univariate L-scale of the marginal uniform random variables. These values further support that default n is large enough. The diagonal is computed from the univariate L-moments of the margins of the Monte Carlo-generated edges and is otherwise uninformative. The secondary diagonal is a rescaling of the \delta^{[\ldots]}_1 by the univariate L-moments of the margins to form L-coscales (nonratios). The copula is perfectly dependent but uncorrelated; so the secondary diagonal has near zeros.

The second L-comoment ratio matrix (coefficient of L-variations and L-correlations) is

  print(theoLClcm$T2)
                [,1]          [,2]
  [1,]  1.000000e+00 -2.079712e-05
  [2,] -2.079712e-05  1.000000e+00

where the diagonal by definition has unities (correlation is unity for a variable on itself) but the secondary diagonal for the L-correlations has near zeros because again the copula is uncorrelated, and the secondary diagonal is computed from the \delta^{[\ldots]}_1. These L-correlations are the Spearman Rho values computed external to the algorithms within rhoCOP.

The third L-comoment ratio matrix (L-skews and L-coskews) is

  print(theoLClcm$T3)
                [,1]          [,2]
  [1,]  3.021969e-06 -2.829783e-05
  [2,] -7.501135e-01  4.518901e-06

where the diagonal by definition should have nero zeros because the univariate L-skew of a uniform variable is zero. These values further support that default n is large enough. The secondary diagonal holds L-coskews. The copula has L-coskew of U wrt V of numerically near zero (symmetry) but measurable asymmetry of L-coskew of V wrt U of \tau^{[21]}_3 \approx -0.75.

The fourth L-comoment ratio matrix (L-kurtosises and L-cokurtosises) is

  print(theoLClcm$T4)
                [,1]          [,2]
  [1,] -2.623665e-06 -3.325177e-05
  [2,] -2.162954e-04 -2.811630e-06

where the diagonal by definition should have nero zeros because the univariate L-kurtosis of a uniform variable is zero—it has no peakedness. These values further support that default n is large enough. The secondary diagonal holds L-cokurtosises and are near zero for this particular copula.

The fifth L-comoment ratio matrix (unnamed) is

  print(theoLClcm$T5)
               [,1]          [,2]
  [1,] 1.813344e-06 -1.296025e-04
  [2,] 1.246436e-01  1.475012e-06

where the diagonal by definition should have nero zeros because the univariate L-kurtosis of a uniform variable is zero—such a random variable has no asymmetry. These values further support that default n is large enough. The secondary diagonal holds the fifth L-comoment ratios. The copula has \tau^{[12]}_5 = 0 of U wrt V of numerically near zero (symmetry) but measurable fifth-order L-comoment asymmetry of V wrt U of \tau^{[21]}_5 \approx 0.125.

Comparison of Sample and Theoretical L-comoments—The previous section shows theoretical values computed as \tau^{[21]}_3 \approx -0.75 and \tau^{[21]}_5 \approx 0.125 for the two L-comoments substantially away from zero. As sample L-comoments these values are samLC$T3[2,1] = \hat\tau^{[21]}_3 \approx -0.751 and samLC$T5[2,1] = \hat\tau^{[21]}_5 \approx 0.123. CONCLUSION: The sample L-comoment algorithms in the lmomco package are validated.

Author(s)

W.H. Asquith

References

Asquith, W.H., 2011, Distributional analysis with L-moment statistics using the R environment for statistical computing: Createspace Independent Publishing Platform, ISBN 978–146350841–8.

Brahimi, B., Chebana, F., and Necir, A., 2015, Copula representation of bivariate L-moments—A new estimation method for multiparameter two-dimensional copula models: Statistics, v. 49, no. 3, pp. 497–521.

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

Serfling, R., and Xiao, P., 2007, A contribution to multivariate L-moments—L-comoment matrices: Journal of Multivariate Analysis, v. 98, pp. 1765–1781.

See Also

lcomCOP, uvlmoms

Examples

## Not run: 
bilmoms(cop=PSP, n=10000, para=NULL, sobol=TRUE)$bilcomoms$T3
# results: Tau3[12]=-0.132, Tau3[21]=-0.132 (Monte Carlo)
lcomCOP(cop=PSP, para=NULL, orders=3)
# results: Tau3[12]=-0.129, Tau3[21]=-0.129 (direct integration)
## End(Not run)

## Not run: 
# This stopped running sometime before June 2023. IS THIS IN COP()?
para <- list(alpha=0.5, beta=0.93, para1=4.5, cop1=GLcop, cop2=PSP)
bilmoms(cop=composite2COP, n=10000, para=para, sobol=TRUE)$bilcomoms$T3
# results: Tau3[12]=0.154, Tau3[21]=-0.0691 (Monte Carlo)
lcomCOP(cop=composite2COP, para=para, orders=3)
# results: Tau3[12]=0.156, Tau3[21]=-0.0668 (direct integration)
## End(Not run)

## Not run: 
UVsim <- simCOP(n=20000, cop=composite2COP, para=para, graphics=FALSE)
samLcom <- lmomco::lcomoms2(UVsim, nmom=5) # sample algorithm
# results: Tau3[12]=0.1489, Tau3[21]=-0.0679 (simulation)
## End(Not run)