R/tests_ArchmChisq.R

In gofCopula: Goodness-of-Fit Tests for Copulae

#' ArchmChisq Gof test using the Anderson-Darling test
#' statistic and the chi-square distribution
#' 
#' \code{\link{gofArchmChisq}} contains the Chisq gof test with a Rosenblatt 
#' transformation for archimedean copulae, described in Hering and Hofert (2015).
#' The test follows the RosenblattChisq test as described in Genest (2009) 
#' and Hofert (2014), and compares the empirical copula against a parametric 
#' estimate of the copula derived under the null hypothesis. The margins can be 
#' estimated by a bunch of distributions and the time which is necessary for 
#' the estimation can be given. The approximate p-values are computed with a 
#' parametric bootstrap, which computation can be accelerated by enabling 
#' in-build parallel computation. The gof statistics are computed with the 
#' function \code{\link{gofTstat}} from the package copula. 
#' It is possible to insert datasets of all dimensions above 1 
#' (except for the \code{"amh"} copula) and the possible copulae are 
#' \code{"clayton"}, \code{"gumbel"}, \code{"frank"}, \code{"joe"} and 
#' \code{"amh"}. The parameter estimation is performed with pseudo maximum
#' likelihood method. In case the estimation fails, inversion of Kendall's tau
#' is used.
#' 
#' This Anderson-Darling test statistic (supposedly) computes
#' U[0,1]-distributed (under \eqn{H_0}{H0}) random variates via the
#' distribution function of chi-square distribution with d degrees of freedom,
#' see Hofert et al. (2014). The \eqn{H_0}{H0} hypothesis is \deqn{C \in
#' \mathcal{C}_0}{C in Ccal0} with \eqn{\mathcal{C}_0}{Ccal0} as the true class
#' of copulae under \eqn{H_0}{H0}.
#' 
#' This test is based on the Rosenblatt transformation for archimedean copula
#' which uses the mapping \eqn{\mathcal{R}: (0,1)^d 
#' \rightarrow (0,1)^d}{R : (0,1)^d -> (0,1)^d}. Following 
#' Hering and Hofert (2015) the mapping provides 
#' pseudo observations \eqn{E_i}{E[i]}, given by \deqn{E_1 = \mathcal{R}(U_1), 
#' \dots, E_n = \mathcal{R}(U_n).}{E_1 = R(U_1), ..., E_n = R(U_n).} 
#' Let \eqn{C} be an Archimedean copula with \eqn{d} monotone generator 
#' \eqn{\psi} and continuous Kendall distribution function \eqn{K_{C}}. Then, 
#' \deqn{e_{j}=\left(\frac{\sum_{k=1}^{j} \psi^{-1}\left(u_{k}
#' \right)}{\sum_{k=1}^{j+1} \psi^{-1}\left(u_{k}\right)}\right)^{j}, 
#' j \in\{1, \ldots, d-1\}}{e[j] = (sum_k^(j)(\psi^(-1)(u_(k))) / 
#' sum_k^(j+1)(\psi^(-1)(u_(k))))^j, j = 1, ..., d-1} and 
#' \deqn{e_{d}= \frac{n}{n+1} K_{C}(C(u))}{e[i] = n / (n + 1) K(C(u))}.
#' 
#' The Anderson-Darling test statistic of the variates
#' 
#' \deqn{G(x_j) = \chi_d^2 \left( x_j \right)}{G(x[j]) = pchisq(x[j], df=d)}
#' 
#' is computed (via \code{ADGofTest::ad.test}), where \eqn{x_j = \sum_{i=1}^d
#' (\Phi^{-1}(e_{ij}))^2}{x[j] =
#' (Phi^{-1}(e_{1j}))^2+...+(Phi^{-1}(e_{dj}))^2}, \eqn{\Phi^{-1}}{Phi^{-1}}
#' denotes the quantile function of the standard normal distribution function,
#' \eqn{\chi_d^2}{pchisq(.,df=d)} denotes the distribution function of the
#' chi-square distribution with \code{d} degrees of freedom, and \eqn{u_{ij}}
#' is the \eqn{j}th component in the \eqn{i}th row of \eqn{\mathbf{u}}{u}.
#' 
#' The test statistic is then given by \deqn{T = -n - \sum_{j=1}^n \frac{2j -
#' 1}{n} [\ln(G(x_j)) + \ln(1 - G(x_{n+1-j}))].}{T = -n - sum((2j - 1)/n
#' [ln(G(x[j])) + ln(1 - G(x[n+1-j]))], j = 1, ..., n).}
#' 
#' The approximate p-value is computed by the formula,
#' 
#' \deqn{\sum_{b=1}^M \mathbf{I}(|T_b| \geq |T|) / M,}{sum(|T[b]| >= |T|, b=1,
#' .., M) / M,}
#' 
#' where \eqn{T} and \eqn{T_b}{T[b]} denote the test statistic and the
#' bootstrapped test statistc, respectively.
#' 
#' For small values of \code{M}, initializing the parallelisation via
#' \code{processes} does not make sense. The registration of the parallel
#' processes increases the computation time. Please consider to enable
#' parallelisation just for high values of \code{M}.
#' 
#' @param copula The copula to test for. Possible are \code{"clayton"}, 
#' \code{"gumbel"}, \code{"frank"}, \code{"joe"} and \code{"amh"}.
#' @param x A matrix containing the data with rows being observations and
#' columns being variables.
#' @param param The copula parameter to use, if it shall not be estimated.
#' @param param.est Shall be either \code{TRUE} or \code{FALSE}. \code{TRUE}
#' means that \code{param} will be estimated.
#' @param margins Specifies which estimation method for the margins shall be
#' used. The default is \code{"ranks"}, which is the standard approach to
#' convert data in such a case. Alternatively the following distributions can
#' be specified: \code{"beta"}, \code{"cauchy"}, Chi-squared (\code{"chisq"}),
#' \code{"f"}, \code{"gamma"}, Log normal (\code{"lnorm"}), Normal
#' (\code{"norm"}), \code{"t"}, \code{"weibull"}, Exponential (\code{"exp"}).
#' Input can be either one method, e.g. \code{"ranks"}, which will be used for
#' estimation of all data sequences. Also an individual method for each margin
#' can be specified, e.g. \code{c("ranks", "norm", "t")} for 3 data sequences.
#' If one does not want to estimate the margins, set it to \code{NULL}.
#' @param flip The control parameter to flip the copula by 90, 180, 270 degrees
#' clockwise. Only applicable for bivariate copula. Default is 0 and possible 
#' inputs are 0, 90, 180, 270 and NULL.
#' @param M Number of bootstrapping loops.
#' @param lower Lower bound for the maximum likelihood estimation of the copula
#' parameter. The constraint is also active in the bootstrapping procedure. The
#' constraint is not active when a switch to inversion of Kendall's tau is
#' necessary. Default \code{NULL}.
#' @param upper Upper bound for the maximum likelihood estimation of the copula
#' parameter. The constraint is also active in the bootstrapping procedure. The
#' constraint is not active when a switch to inversion of Kendall's tau is
#' necessary. Default \code{NULL}.
#' @param seed.active Has to be either an integer or a vector of M+1 integers.
#' If an integer, then the seeds for the bootstrapping procedure will be
#' simulated. If M+1 seeds are provided, then these seeds are used in the
#' bootstrapping procedure. Defaults to \code{NULL}, then \code{R} generates
#' the seeds from the computer runtime. Controlling the seeds is useful for
#' reproducibility of a simulation study to compare the power of the tests or
#' for reproducibility of an empirical study.
#' @param processes The number of parallel processes which are performed to
#' speed up the bootstrapping. Shouldn't be higher than the number of logical
#' processors. Please see the details.
#' @return An object of the \code{class} gofCOP with the components
#' \item{method}{a character which informs about the performed analysis}
#' \item{copula}{the copula tested for} \item{margins}{the method used to
#' estimate the margin distribution.} \item{param.margins}{the parameters of
#' the estimated margin distributions. Only applicable if the margins were not
#' specified as \code{"ranks"} or \code{NULL}.} \item{theta}{dependence
#' parameters of the copulae} \item{df}{the degrees of freedem of the copula.
#' Only applicable for t-copula.} \item{res.tests}{a matrix with the p-values
#' and test statistics of the hybrid and the individual tests}
#' @references Christian Genest, Bruno Remillard, David Beaudoin (2009).
#' Goodness-of-fit tests for copulas: A review and a power study.
#' \emph{Insurance: Mathematics and Economics, Volume 44, Issue 2, April 2009,
#' Pages 199-213, ISSN 0167-6687}.
#' \doi{10.1016/j.insmatheco.2007.10.005}\cr \cr Marius
#' Hofert, Ivan Kojadinovic, Martin Maechler, Jun Yan (2014). copula:
#' Multivariate Dependence with Copulas. \emph{R package version 0.999-15.}.
#' \url{https://cran.r-project.org/package=copula} \cr \cr
#' Christian Hering, Marius Hofert (2015) Goodness-of-fit tests for 
#' Archimedean copulas in high dimensions. \emph{In: Glau K., Scherer M., 
#' Zagst R. (eds) Innovations in Quantitative Risk Management, 
#' Springer Proceedings in Mathematics & Statistics, Volume 99, 
#' Springer, Cham, 357-373}.
#' \doi{10.1007/978-3-319-09114-3_21} \cr \cr
#' @examples
#' 
#' data(IndexReturns2D)
#' 
#' gofArchmChisq("clayton", IndexReturns2D, M = 10)
#' 
#' @export gofArchmChisq
gofArchmChisq <- function(copula = c("clayton", "gumbel", "frank", "joe", 
                                     "amh"), x, param = 0.5, param.est = TRUE, 
                          margins = "ranks", flip = 0, M = 1000, lower = NULL, 
                          upper = NULL, seed.active = NULL, processes = 1) {
  if (is.matrix(x) == FALSE) {
    stop("x must be a matrix")
  }
  if (length(copula) > 1) {
stop(
"'copula' has to be a vector of length 1. Please select only one copula."
)
  }
  if (is.element(copula, c("clayton", "gumbel", "frank", "joe", 
                           "amh")) == FALSE) {
    stop("This copula is not implemented for gofArchmChisq.")
  }
  if (!is.numeric(processes)) {
    stop("The argument 'processes' has to be a numeric.")
  }
  if (processes %% 1 != 0 | processes < 1) {
    stop("The argument 'processes' has to be a positive integer.")
  }
  if (!is.numeric(M)) {
    stop("The argument 'M' has to be a numeric.")
  }
  if (M %% 1 != 0 | M < 0) {
    stop("The argument 'M' has to be a positive integer.")
  }
  if (!is.numeric(param)) {
    stop("The argument 'param' has to be a numeric.")
  }
  if (!inherits(param.est, "logical")) {
    stop("The argument 'param.est' has to be either 'TRUE' or 'FALSE'.")
  }
  if (!is.null(lower) & !is.numeric(lower) | !is.null(upper) & 
      !is.numeric(upper)) {
    stop("The arguments 'upper' and 'lower' have to be either NULL or numeric.")
  }
  if (copula == "gumbel" & param.est == FALSE & param <= 1) {
    param <- 1.5
warning(
"When copula is 'gumbel', 'param' has to be larger 1. Because 'param.est' was 
set to 'FALSE', 'param' was set to 1.5 as default value for 'gumbel' copula."
)
  }
  if (!is.null(seed.active) & length(seed.active) != 1 & 
      length(seed.active) != (M + 1)) {
    stop("The seed has to be an integer or a vector of M+1 seeds.")
  }
  if (!is.null(seed.active) & length(seed.active) == 1) {
    set.seed(seed.active)
    RNGsetting <- RNGkind()
    RNGkind(sample.kind = "default")
    on.exit(RNGkind(sample.kind = RNGsetting[3]))
    seed.active <- sample(x = 2147483647, size = M + 1)
  }
  if (!is.null(seed.active) & all(!vapply(seed.active, 
                                          function(x) x %% 1 == 0, TRUE))) {
stop(
"All seeds have to be whole numbers. Please check seed.active for non-whole 
numbers."
)
  }
  
  # estimation of margins and copula parameters
  erg <- .margins.param.est(copula = copula, margins = margins, x = x, 
                            param = param, param.est = param.est, df = 4, 
                            df.est = TRUE, dispstr = "ex", lower = lower, 
                            upper = upper, flip = flip)
  # test with parametric bootstrap. Switch to Kendall's Tau if Maximum
  # Likelihood estimation fails
  res <- try(.gofCopulapb(copula = erg[[1]], x = erg[[2]], M = M, 
                          method = "ArchmAnChisq", estim.method = "mpl", 
                          processes = processes, param.est = param.est, 
                          df.est = erg[[5]], dispstr = "ex", 
                          param.margins = erg[[4]], margins = margins, 
                          seed.active = seed.active, lower = lower, 
                          upper = upper, flip = erg[[6]]), silent = TRUE)
  if (inherits(res, "try-error")) {
warning(
"Pseudo Maximum Likelihood estimation of the parameters while the bootstrapping 
procedure failed. The estimation was performed with inversion of Kendall's Tau. 
Therefore df.est was set to FALSE for the bootstrapping."
)
    res.f <- .gofCopulapb(copula = erg[[1]], x = erg[[2]], M = M, 
                          method = "ArchmAnChisq", estim.method = "itau", 
                          processes = processes, param.est = param.est, 
                          df.est = FALSE, dispstr = "ex", 
                          param.margins = erg[[4]], margins = margins, 
                          seed.active = seed.active, lower = lower, 
                          upper = upper, flip = erg[[6]])
    return(res.f)
  } else {
    return(res)
  }
}