# Bivariate L-moments and L-comoments of a Copula 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 [REMOVE_ME]$$ \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{,} [REMOVE_ME_2]$$ [REMOVE_ME]$$ \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{,} [REMOVE_ME_2]$$ [REMOVE_ME]$$ \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} [REMOVE_ME_2]$$ [REMOVE_ME]$$ \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{,} [REMOVE_ME_2]$$ where the bivariate L-moments are related to the L-comoment ratios by [REMOVE_ME]$$ 6\delta^{[12]}_k = \tau^{[12]}_{k+1}\mbox{\quad and \quad}6\delta^{[21]}_k = \tau^{[21]}_{k+1}\mbox{,} [REMOVE_ME_2]$$ where in otherwords, the third bivariate L-moment delta^{[12]}_3 is one sixth theL-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$. utf8 ## 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 theL-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$. ```r 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. ## Returns An `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 `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 `NA`s. 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.** ## 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. ## Author(s) W.H. Asquith ## See Also `lcomCOP`, `uvlmoms` ## Examples ```r ## 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) ```

bilmoms function