copula (version 1.1-3)

RSpobs: Pseudo-Observations of Radial and Uniform Part of Elliptical and Archimedean Copulas

Description

Given a matrix of iid multivariate data from a meta-elliptical or meta-Archimedean model, RSpobs() computes pseudo-observations of the radial part \(R\) and the vector \(\bm{S}\) which follows a uniform distribution on the unit sphere (for elliptical copulas) or the unit simplex (for Archimedean copulas). These quantities can be used for (graphical) goodness-of-fit tests, for example.

Usage

RSpobs(x, do.pobs = TRUE, method = c("ellip", "archm"), ...)

Value

A list with components R (an

\(n\)-vector containing the pseudo-observations of the radial part) and S (an \((n, d)\)-matrix containing the pseudo-observations of the uniform distribution (on the unit sphere/simplex)).

Arguments

x

an \((n, d)\)-matrix of data; if do.pobs=FALSE, the rows of x are assumed to lie in the \(d\)-dimensional unit hypercube (if they do not, this leads to an error).

do.pobs

logical indicating whether pobs() is applied to x for transforming the data to the \(d\)-dimensional unit hypercube.

method

character string indicating the assumed underlying model, being meta-elliptical if method="ellip" (in which case S should be approximately uniform on the \(d\)-dimensional unit sphere) or meta-Archimedean if method="archm" (in which case S should be approximately uniform on the \(d\)-dimensional unit simplex).

...

additional arguments passed to the implemented methods. These can be

method="ellip"

qQg() (the quantile function of the (assumed) distribution function \(G_g\) as given in Genest, Hofert, G. Nešlehová (2014)); if provided, qQg() is used in the transformation for obtaining pseudo-observations of \(R\) and \(\bm{S}\) (see the code for more details).

method="archm"

iPsi() (the assumed underlying generator inverse); if provided, iPsi() is used in the transformation for obtaining pseudo-observations of \(R\) and \(\bm{S}\) (see the code for more details).

Details

The construction of the pseudo-obersvations of the radial part and the uniform distribution on the unit sphere/simplex is described in Genest, Hofert, G. Nešlehová (2014).

References

Genest, C., Hofert, M., G. Nešlehová, J., (2014). Is the dependence Archimedean, elliptical, or what? To be submitted.

See Also

pobs() for computing the “classical” pseudo-observations.

Examples

set.seed(100)
n <- 250 # sample size
d <- 5 # dimension
nu <- 3 # degrees of freedom

## Build a mean vector and a dispersion matrix,
## and generate multivariate t_nu data:
mu <- rev(seq_len(d)) # d, d-1, .., 1
L <- diag(d) # identity in dim d
L[lower.tri(L)] <- 1:(d*(d-1)/2)/d # Cholesky factor (diagonal > 0)
Sigma <- crossprod(L) # pos.-def. dispersion matrix (*not* covariance of X)
X <- rep(mu, each=n) + mvtnorm::rmvt(n, sigma=Sigma, df=nu) # multiv. t_nu data
## note: this is *wrong*: mvtnorm::rmvt(n, mean=mu, sigma=Sigma, df=nu)

## compute pseudo-observations of the radial part and uniform distribution
## once for 1a), once for 1b) below
RS.t    <- RSpobs(X, method="ellip", qQg=function(p) qt(p, df=nu)) # 'correct'
RS.norm <- RSpobs(X, method="ellip", qQg=qnorm) # for testing 'wrong' distribution
stopifnot(length(RS.norm$R) == n, length(RS.t$R) == n,
          dim(RS.norm$S) == c(n,d), dim(RS.t$S) == c(n,d))

## 1) Graphically testing the radial part

## 1a) Q-Q plot of R against the correct quantiles
qqplot2(RS.t$R, qF=function(p) sqrt(d*qf(p, df1=d, df2=nu)),
        main.args=list(text= substitute(bold(italic(F[list(d.,nu.)](r^2/d.))~~"Q-Q Plot"),
                                        list(d.=d, nu.=nu)),
		       side=3, cex=1.3, line=1.1, xpd=NA))

## 1b) Q-Q plot of R against the quantiles of F_R for a multivariate normal
##     distribution
qqplot2(RS.norm$R, qF=function(p) sqrt(qchisq(p, df=d)),
        main.args=list(text= substitute(bold(italic(chi[D_]) ~~ "Q-Q Plot"), list(D_=d)),
               side=3, cex=1.3, line=1.1, xpd=NA))

## 2) Graphically testing the angular distribution

## auxiliary function
qqp <- function(k, Bmat) {
    d <- ncol(Bmat) + 1
    qqplot2(Bmat[,k],
            qF = function(p) qbeta(p, shape1=k/2, shape2=(d-k)/2),
            main.args=list(text= substitute(plain("Beta")(s1,s2) ~~ bold("Q-Q Plot"),
                                            list(s1 = k/2, s2 = (d-k)/2)),
  	                   side=3, cex=1.3, line=1.1, xpd=NA))
}

## 2a) Q-Q plot of the 'correct' angular distribution
##     (Bmat[,k] should follow a Beta(k/2, (d-k)/2) distribution)
Bmat.t <- gofBTstat(RS.t$S)
qqp(1, Bmat=Bmat.t) # k=1
qqp(3, Bmat=Bmat.t) # k=3

## 2b) Q-Q plot of the 'wrong' angular distribution
Bmat.norm <- gofBTstat(RS.norm$S)
qqp(1, Bmat=Bmat.norm) # k=1
qqp(3, Bmat=Bmat.norm) # k=3

## 3) Graphically check independence between radial part and B_1 and B_3

## 'correct' distributions (multivariate t)
plot(pobs(cbind(RS.t$R, Bmat.t[,1])), # k = 1
          xlab=quote(italic(R)), ylab=quote(italic(B)[1]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[1])))
plot(pobs(cbind(RS.t$R, Bmat.t[,3])), # k = 3
	  xlab=quote(italic(R)), ylab=quote(italic(B)[3]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[3])))

## 'wrong' distributions (multivariate normal)
plot(pobs(cbind(RS.norm$R, Bmat.norm[,1])), # k = 1
          xlab=quote(italic(R)), ylab=quote(italic(B)[1]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[1])))
plot(pobs(cbind(RS.norm$R, Bmat.norm[,3])), # k = 3
	  xlab=quote(italic(R)), ylab=quote(italic(B)[3]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[3])))