Attention: This function is thought to be mostly 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, ...)A copula function;
Vector of parameters or other data structure, if needed, to pass to the copula;
A logical to trigger return of the \(\delta_k\) and skip L-comoment computation;
The Monte Carlo integration size. The default seems to be at least an order of magnitude greater than needed for many applied problems;
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; and
Additional arguments to pass to the densityCOP function.
An R list of the bivariate L-moments is returned.
The bivariate L-moments \(\delta^{[12]}_k\) of \(U\) with respect to \(V\) for \(k \in [1,2,3,4]\);
The bivariate L-moments \(\delta^{[21]}_k\) of \(V\) with respect to \(U\) for \(k \in [1,2,3,4]\);
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;
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
An attribute identifying the computational source of the bivariate L-moments and bivariate L-comoments: “bilmoms.”
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.
# NOT RUN {
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.155, Tau3[21]=-0.0667 (Monte Carlo)
lcomCOP(cop=composite2COP, para=para, orders=3)
# results: Tau3[12]=0.156, Tau3[21]=-0.0668 (direct integration)
# }
# 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)
# }
Run the code above in your browser using DataLab