R/07_choose_copula.R

In rdmatheus/BCNSM: Multivariate Box-Cox Symmetric Distributions Generated by a Normal Scale Mixture Copula

#' Estimation of the Extra Parameter of the NSM Copula in the Fit of the BCNSM Distribution
#'
#' Estimation of the extra parameter of the NSM copula in the fit of the BCNSM distribution
#' via profile log-likelihood.
#'
#' @param object an \code{"bcnsm"} object.
#' @param grid grid of values that will be used to evaluate the profiled log-likelihood function.
#'     If the tweaks are not computationally intensive, we suggest \code{grid = 1:15}.
#' @param copula character; informs which normal scale mixture distribution
#'     should be used to generate the NSM copula. Currently,
#'     the copulas available are: Gaussian (\code{"gaussian"}),
#'     Student's t (\code{"t"}), slash (\code{"slash"}), and hyperbolic (\code{"hyp"}).
#' @param trace logical; if \code{TRUE}, a summary with the profiled log-likelihood value, the AIC,
#'     the BIC, and the run-time of the fit is displayed.
#' @param plot logical; if \code{TRUE}, a graph of the profiled log-likelihood evaluated in the
#'     considered grid of values is shown.
#' @param control a list of control arguments specified via \code{\link[bcnsm]{control_fit}}.
#' @param ... further arguments passed to \code{\link[bcnsm]{control_fit}}.
#'
#' @return An object of class \code{"choose_copula"}. More specifically, it returns a list in which
#'     each element consists of the fit of the BCNSM distribution with each value of the extra
#'     parameter specified in \code{grid}.
#'
#' @references Vanegas, L. H., and Paula, G. A. (2016). Log-symmetric distributions: statistical properties and
#'     parameter estimation. *Brazilian Journal of Probability and Statistics*, 30, 196-220.
#'
#' Ferrari, S. L. P., and Fumes, G. (2017). Box-Cox symmetric distributions and applications to
#'     nutritional data. *AStA Advances in Statistical Analysis*, 101, 321-344.
#'
#' Medeiros, R. M. R. de, and Ferrari, S. L. P. (2023). Multivariate Box-Cox symmetric distributions
#'     generated by a normal scale mixture copula.
#'
#' @author Rodrigo M. R. de Medeiros <\email{rodrigo.matheus@live.com}>
#'
#' @export
#'
choose_copula <- function(object, grid, copula, trace = TRUE, plot = TRUE,
                          control = control_fit(...), ...) {
  n <- object$nobs
  fit_update <- lapply(grid, function(delta) {
    init <- Sys.time()
    opt <- try(stats::update(object, copula = copula, delta = delta))
    end <- Sys.time()

    if (trace) {
      cat(
        "\ndelta:", delta,
        "|",
        "logLik:", if (inherits(opt, "error")) NA else round(stats::logLik(opt), 3),
        "|",
        "AIC:", if (inherits(opt, "error")) NA else round(stats::AIC(opt), 3),
        "|",
        "BIC:", if (inherits(opt, "error")) NA else round(stats::AIC(opt, k = log(n)), 3),
        "|",
        "Time to run:", if (inherits(opt, "error")) NA else round(difftime(end, init, units = "secs"), 3), "secs"
      )
    }

    opt
  })

  if (plot) {

    ll <- vector()
    for (i in 1:length(grid)) {
      ll[i] <- if (inherits(fit_update[[i]], "error")) NA else stats::logLik(fit_update[[i]])
    }

    plot(grid, ll, type = "o", pch = 16, cex = 0.6,
         xlab = expression(delta), ylab = "Profile log-likelihood")
    graphics::abline(v = grid[which.max(ll)], lty = 3, col = "grey", lwd = 2)
    graphics::points(c(grid[which.max(ll)], grid[which.max(ll)]), c(ll[which.max(ll)], ll[which.max(ll)]),
                     col = c("#56B1F7", 1), pch = c(16, 1))


  }

  fit_update$copula <- copula
  fit_update$grid <- grid
  class(fit_update) <- "choose_copula"
  fit_update
}

#' @name choose_copula-methods
#' @title Methods for 'choose_copula' objects
#' @param x an object of class \code{"choose_copula"}.
#' @param ... further arguments passed to or from other methods.
#'
#' @author Rodrigo M. R. de Medeiros <\email{rodrigo.matheus@live.com}>
NULL

# Print
#' @rdname choose_copula-methods
#' @export
print.choose_copula <- function(x, ...) {

  cat(crayon::cyan("BCNSM fit with", x$copula, "copula\n"))

  grid <- x$grid
  n <- x[[1]]$nobs
  i <- 1
  ll <- AIC <- BIC <- vector("numeric", length(grid))
  for (delta in grid) {
    ll[i] <- if (inherits(x[[i]], "error")) NA else as.numeric(stats::logLik(x[[i]]))
    AIC[i] <- if (inherits(x[[i]], "error")) NA else round(stats::AIC(x[[i]]), 3)
    BIC[i] <- if (inherits(x[[i]], "error")) NA else round(stats::AIC(x[[i]], k = log(n)), 3)

    cat(
      "\ndelta:", delta,
      "|",
      "logLik:", round(ll[i], 3),
      "|",
      "AIC:", AIC[i],
      "|",
      "BIC:", BIC[i]
    )

    i <- i + 1
  }

  cat(crayon::cyan("\n\nBest value for delta according to AIC:"), grid[which.min(AIC)], crayon::cyan("BIC:"), grid[which.min(BIC)],
      crayon::cyan("and logLik:"), grid[which.max(ll)])

  invisible(x)
}


# Plot
#' @rdname choose_copula-methods
#' @export
plot.choose_copula <- function(x, ...) {
  grid <- x$grid

  ll <- vector()
  for (i in 1:length(grid)) {
    ll[i] <- if (inherits(x[[i]], "error")) NA else stats::logLik(x[[i]])
  }

  plot(grid, ll, type = "o", pch = 16, cex = 0.6,
       xlab = expression(delta), ylab = "Profile log-likelihood")
  graphics::abline(v = grid[which.max(ll)], lty = 3, col = "grey", lwd = 2)
  graphics::points(c(grid[which.max(ll)], grid[which.max(ll)]), c(ll[which.max(ll)], ll[which.max(ll)]),
                   col = c("#56B1F7", 1), pch = c(16, 1))

  invisible(x)
}