| bilmoms | R Documentation |
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.
bilmoms(cop=NULL, para=NULL, only.bilmoms=FALSE, n=1E5,
sobol=TRUE, scrambling=0, ...)
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 |
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 |
scrambling |
The argument of the same name for |
... |
Additional arguments to pass to the |
An R list of the bivariate L-moments is returned.
bilmomUV |
The bivariate L-moments |
bilmomVU |
The bivariate L-moments |
error.rho |
An “error” term in units of |
bilcomoms |
If not |
source |
An attribute identifying the computational source of the bivariate L-moments and bivariate L-comoments: “bilmoms.” |
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.
W.H. Asquith
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.
lcomCOP, uvlmoms
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.