RAYcop {copBasic}

The Rayleigh Copula

Description

The Rayleigh copula (Boškoskia and others, 2018) is

\mathbf{C}_{\Theta}(u,v) = \mathbf{RAY}(u,v; \Theta) = 1 + A - B\mbox{,}

A = e^{\Theta a_2 - a_2}\biggl(e^{-a_1}\int_0^{\Theta a_2} e^{-s}I_0\bigl(2\sqrt{a_1 s}\bigr)\,\mathrm{d}s - 1\biggr)\,\mathrm{d}s\mbox{,}

B = e^{-a_1}\int_0^{a_2}e^{-s}I_0\bigl(2\sqrt{\Theta a_1 s}\bigr)\,\mathrm{d}s\mbox{,}

where a1 = -\log(1-u)/(1-\Theta), a2 = -\log(1-v)/(1-\Theta), I_\nu(x) is the modified Bessel function of the first kind of order \nu (see base::besselI()), and \Theta \in (0,1]. The copula, as \Theta \rightarrow 0^{+} limits, to the independence coupla (\mathbf{\Pi}(u,v); P) and as \Theta \rightarrow 1^{-} limits to the comonotonicity copula (\mathbf{M}(u,v); M). Finally, there are formulations of the Rayleigh copula using the Marcum-Q function, but the copBasic developer has not been able to make such work. If the Marcum-Q function could be used, then only one integration and not the two involving the modified Bessel function are possible. Infinite integrations begin occurring in the upper right corner for about \Theta > 0.995 at which point the \mathbf{M}(u,v) copula is called in the source code.

Usage

RAYcop(u, v, para=NULL, rho=NULL, method=c("default"),
             rel.tol=.Machine$double.eps^0.5, ...)

Arguments

u

Nonexceedance probability u in the X direction;

v

Nonexceedance probability v in the Y direction;

para

A vector (single element) of parameters—the \Theta parameter of the copula;

rho

Value for Spearman Rho from which parameter \Theta is computed by polynomial approximation and returned. The estimation appears sufficient for most pratical applications (see Examples);

method

The computational method of integrals associated with the definition of the copula; this is designed for the ability to switch eventually in sources to Marcum-Q function implementation. The definition in January 2023 and default is to call the two Bessel function integrals shown for the definition in this documentation;

rel.tol

Argument of the same name for integrate() call; and

...

Additional arguments to pass.

Value

Value(s) for the copula are returned.

Note

The documentation in Zeng and other [Part II] appear to have corrected the Marcum-Q function solution to the copula. The essence of that solution is with a Chi-distribution computation the Marcum-Q. Testing indictates that this is a correct solution, but the derivative for the conditional simulation as built into the design of copBasic has difficulties. Perhaps this is related to numerical precision of the Marcum-Q?

  sapply(seq_len(length((u))), function(i) {
    a1 <- -log(1-u[i])
    if(is.infinite(a1)) return(v[i])
    a2 <- -log(1-v[i])
    if(is.infinite(a2)) return(u[i])
    a1 <- exp(log(a1) - log(1-p))
    a2 <- exp(log(a2) - log(1-p))
    a3 <- marcumq.chi(sqrt(2*  a1), sqrt(2*p*a2)) # Zeng and others (Part II)
    a4 <- marcumq.chi(sqrt(2*p*a1), sqrt(2*  a2)) # Zeng and others (Part II)
    zz <- 1 + (1-v[i])*a3 - (1-u[i])*(1-a4)       # Zeng and others (Part II)
    zz[zz < 0] <- 0
    zz[zz > 1] <- 1
    return(zz)
  })

Author(s)

W.H. Asquith

References

Boškoskia, P., Debenjaka, A., Boshkoskab, B.M., 2018, Rayleigh copula for describing impedance data with application to condition monitoring of proton exchange membrane fuel cells: European Journal of Operational Research, v. 266, pp. 269–277, doi:10.1016/j.ejor.2017.08.058.

Zeng, X., Ren, J., Wang, Z., Marshall, S., and Durrani, T., [undated], Copulas for statistical signal processing (Part I)—Extensions and generalization, accessed January 14, 2024, at https://pure.strath.ac.uk/ws/portalfiles/portal/34078849/Copulas_Part1_v2_6.pdf.

Zeng, X., Ren, J., Sun, M., Marshall, S., and Durrani, T., [undated], Copulas for statistical signal processing (Part II)—Simulation, optimal selection and practical applications, accessed January 14, 2024, at https://strathprints.strath.ac.uk/48371/1/Copulas_Part2s_v2_5_2.pdf

See Also

M, P

Examples

RAYcop(0.2, 0.8, para=0.8) # [1] 0.1994658  (by the dual Bessel functions)

RAYcop(0.8, 0.2, para=RAYcop(rho=rhoCOP(cop=RAYcop, para=0.8)))
# [1] 0.1994651 from polynomial conversion of Rho to Theta

## Not run: 
# Recipe for assembling the Spearman Rho to Theta polynomial in sources.
Thetas <- seq(0, 0.999, by=0.001); RHOs <- NULL
for(p in Thetas) RHOs <- c(RHOs, rhoCOP(cop=RAYcop, para=p))
LM <- lm(Thetas ~ RHOs + I(RHOs^2) + I(RHOs^4) + I(RHOs^6) - 1 )
Rho2Theta <- function(rho) {
  coes <- c(1.32682824, -0.38876290, 0.09072305, -0.02921836)
  sapply(rho, function(r) coes[1]*r^1 + coes[2]*r^2 + coes[3]*r^4 + coes[4]*r^6 )
}
plot(RHOs, Thetas, type="l", col=grey(0.8), lwd=12, lend=1,
      xlab="Spearman Rho", ylab="Rayleigh Copula Parameter Theta")
lines(RHOs, Rho2Theta(RHOs), col="red", lwd=2) # 
## End(Not run)