Compute the Kullback--Leibler Divergence, Jeffrey Divergence, and Kullback--Leibler sample size following Joe (2014, pp. 234--237). Consider two densities \(f = c_1(u,v; \Theta_f)\) and \(g = c_2(u,v; \Theta_g)\) for two different bivariate copulas \(\mathbf{C}_1(\Theta_1)\) and \(\mathbf{C}_2(\Theta_2)\) having respective parameters \(\Theta\), then the Kullback--Leibler Divergence of \(f\) relative to \(g\) is $$\mathrm{KL}(f {\mid} g) = \int\!\!\int_{\mathcal{I}^2} g\, \log(g/f)\,\mathrm{d}u\mathrm{d}v\mbox{,}$$ and Kullback--Leibler Divergence of \(g\) relative to \(f\) is $$\mathrm{KL}(g {\mid} f) = \int\!\!\int_{\mathcal{I}^2} f\, \log(f/g)\,\mathrm{d}u\mathrm{d}v\mbox{,}$$ where the limits of integration \(\mathcal{I}^2\) theoretically are closed on \([0,1]^2\) but an open interval \((0,1)^2\) might be needed for numerical integration. Note, in general \(\mathrm{KL}(f {\mid} g) \ne \mathrm{KL}(g {\mid} f)\). The \(\mathrm{KL}(f {\mid} g)\) is the expected log-likelihood ratios of \(g\) to \(f\) when \(g\) is the true density (Joe, 2014, p. 234), whereas \(\mathrm{KL}(g {\mid} f)\) is the opposite.
This asymmetry leads to Jeffrey Divergence, which is defined as a symmetrized version of the two Kullback--Leibler Divergences, and is $$J(f,g) = \mathrm{KL}(f {\mid} g) + \mathrm{KL}(g {\mid} f) = \int\!\!\int_{\mathcal{I}^2} (g-f)\, \log(g/f)\,\mathrm{d}u\mathrm{d}v\mbox{.}$$
The variances of the Kullback--Leibler Divergences are defined as $$\sigma^2_{\mathrm{KL}(f {\mid} g)} = \int\!\!\int_{\mathcal{I}^2} g\,[\log(g/f)]^2\,\mathrm{d}u\mathrm{d}v - [\mathrm{KL}(f|g)]^2\mbox{,}$$ and $$\sigma^2_{\mathrm{KL}(g {\mid} f)} = \int\!\!\int_{\mathcal{I}^2} f\,[\log(f/g)]^2\,\mathrm{d}u\mathrm{d}v - [\mathrm{KL}(g|f)]^2\mbox{.}$$
For comparison of copula families \(f\) and \(g\) and taking an \(\alpha = 0.05\), the Kullback--Leibler sample size is defined as $$n_{f\!g} = \bigl[\Phi^{(-1)}(1-\alpha) \times \eta_\mathrm{KL}\bigr]^2\mbox{,}$$ where \(\Phi^{(-1)}(t)\) is the quantile function for the standard normal distribution \(\sim\) N(0,1) for nonexceedance probability \(t\), and \(\eta_\mathrm{KL}\) is the maximum of $$\eta_\mathrm{KL} = \mathrm{max}\bigl[\sigma_{\mathrm{KL}(f {\mid} g)}/\mathrm{KL}(f {\mid} g),\, \sigma_{\mathrm{KL}(g {\mid} f)}/\mathrm{KL}(g {\mid} f)\bigr]\mbox{.}$$ The \(n_{f\!g}\) gives an indication of the sample size needed to distinguish \(f\) and \(g\) with a probability of at least \(1 - \alpha = 1 - 0.05 = 0.95\) or 95 percent.
The copBasic features a naïve Monte Carlo integration scheme in the primary interface kullCOP, although the function kullCOPint provides for nested numerical integration. This later function is generally fast but suffers too much for general application from integral divergencies issued from the integrate() function in R---this must be judged in the light that the copBasic package focuses only on implementation of the function of the copula itself and numerical estimation of copula density (densityCOP) and not analytical copula densities or hybrid representations thereof. Sufficient “bread crumbs” are left among the code and documentation for users to re-implement if speed is paramount. Numerical comparison to the results of Joe (2014) (see Examples) suggests that the default Monte Carlo sample size should be more than sufficient for general inference with the expense of considerable CPU time; however, a couple of repeated calls of kullCOP would be advised and compute say the mean of the resulting sample sizes.
kullCOP(cop1=NULL, cop2=NULL, para1=NULL, para2=NULL, alpha=0.05,
del=0, n=1E5, verbose=TRUE, sobol=FALSE, scrambling=0, ...)kullCOPint(cop1=NULL, cop2=NULL, para1=NULL, para2=NULL, alpha=0.05,
del=.Machine$double.eps^0.25, verbose=TRUE, ...)
An R
list is returned having the following components:
kullCOP (Monte Carlo integration) only---The simulation size for numerical integration;
A vector of the Kullback--Leibler Divergences and their standard deviations: \(\mathrm{KL}(f {\mid} g)\), \(\sigma_{\mathrm{KL}(f {\mid} g)}\), \(\mathrm{KL}(g {\mid} f)\), and \(\sigma_{\mathrm{KL}(g {\mid} f)}\), respectively;
kullCOP (Monte Carlo integration) only---The standard deviation of the divergences and the variances;
Jeffrey Divergence \(J(f,g)\);
Kullback--Leibler sample size \(n_{f\!g}\); and
kullCOPint (numerical integration) only---An R list of the outer call of the integrate() function for the respective numerical integrals shown in this documentation.
A copula function corresponding to copula \(f\) in Joe (2014);
Vector of parameters or other data structure, if needed, to pass to the copula \(f\);
A copula function corresponding to copula \(g\) in Joe (2014);
Vector of parameters or other data structure, if needed, to pass to the copula \(g\);
The \(\alpha\) in the Kullback--Leibler sample size equation;
A small value used to denote the lo and hi values of the numerical integration: lo = del and hi = 1 - del. If del == 0, then lo = 0 and hi = 1, which corresponds to the theoretical limits \(\mathcal{I}^2 = [0,1]^2\) and are defaulted here to \([0,1]^2\) because the Monte Carlo algorithm is preferred for general application. The end point control, however, is maintained just in case pathological situations should arise;
kullCOP (Monte Carlo integration) only---the Monte Carlo integration simulation size;
A logical trigging a couple of status lines of output through the message() function in R;
A logical trigging Sobol sequences for the Monte Carlo integration instead of the bivariate uniform distribution. The Sobol sequences are dependent on the randtoolbox package and 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, the randtoolbox at least at version 2.0.+ has “scrambling” of Sobol sequences temporarily disabled, and hence scrambling=0 as default for kullCOP;
The argument of the same name for randtoolbox::sobol; and
Additional arguments to pass to the densityCOP function.
W.H. Asquith
Joe, H., 2014, Dependence modeling with copulas: Boca Raton, CRC Press, 462 p.
densityCOP, vuongCOP