R/estimation.R

In rdmatheus/mbcsec: Multivariate Box-Cox Symmetric Distributions Generated From Elliptical Copula

## Options for optim function --------------------------------------------------
control_fit <- function(method = "BFGS", maxit = 100, hessian = FALSE,
                    trace = FALSE, flambda = FALSE, m_start = NULL,
                    P_start = NULL, fcopula = FALSE, c_start = NULL, ...)
{
  rval <- list(flambda = flambda, method = method, maxit = maxit,
               hessian = hessian, trace = trace, m_start = m_start,
               P_start = P_start, fcopula = fcopula, c_start = c_start)

  rval <- c(rval, list(...))

  if (!is.null(rval$fnscale))
    warning("fnscale must not be modified")

  rval$fnscale <- -1

  if (is.null(rval$reltol))
    rval$reltol <- .Machine$double.eps^(1/1.2)

  rval
}

## Marginal Initial values -----------------------------------------------------
m_inits <- function(data, link = list(link.mu = "log", link.sigma = "log"),
                    flambda = FALSE, gen = "NO"){

  y <- data$y
  X <- data$X
  Z <- data$Z
  S <- data$S
  n <- length(y)

  if (is.null(X)) X <- matrix(1, nrow = n)
  if (is.null(Z)) Z <- matrix(1, nrow = n)
  if (is.null(S)) S <- matrix(1, nrow = n)

  k1 <- ncol(X); k2 <- ncol(Z); k3 <- ncol(S)

  link.mu <- stats::make.link(link$link.mu)
  link.sigma <- stats::make.link(link$link.sigma)

  ### Marginal initial values
  beta0 <- as.numeric(solve(t(X)%*%X)%*%t(X)%*%link.mu$linkfun(y))
  gamma0 <- c(link.sigma$linkfun(stats::sd(y)/mean(y)), rep(0, NCOL(Z) - 1))
  zeta0 <- nu0 <- names_zeta <- names_nu <- NULL

  if (!flambda){
    zeta0 <- c(-e1071::skewness(y), rep(0, NCOL(S) - 1))
    names_zeta <- paste("zeta", 1:k3, sep = "")
  }

  if (BCSgen(gen)$extraP == TRUE){
    nu0 <- e1071::kurtosis(y) + 3
    names_nu <- "nu"
  }

  m_start <- c(beta0, gamma0, zeta0, nu0)
  names(m_start) <- c(paste("beta", 1:k1, sep = ""),
                      paste("gamma", 1:k2, sep = ""),
                      names_zeta,
                      names_nu)
  m_start
}

## Direct maximization for marginal regression ---------------------------------
mreg_fit <- function(data, cormat = ind(),
                     link = list(link.mu = "log", link.sigma = "log"),
                     copula = "gaussian", gen = "NO",
                     control = control_fit(...), ...){

  ### Data setting -------------------------------------------------------------
  y <- data$y
  X <- data$X
  Z <- data$Z
  S <- data$S
  n <- length(y)

  ### Control list -------------------------------------------------------------
  flambda <- control$flambda
  method <- control$method
  maxit <- control$maxit
  hessian <- control$hessian
  trace <- control$trace
  m_start  <- control$m_start
  P_start <- control$P_start
  fcopula <- control$fcopula
  c_start <- control$c_start
  control$flambda <- control$method <- control$hessian <- control$m_start <-
  control$P_start <- control$c_start <- control$fcopula <- NULL

  ### Initial values -----------------------------------------------------------

  #### Margins -----------------------------------------------------------------
  if (is.null(X)) X <- matrix(1, nrow = n)
  if (is.null(Z)) Z <- matrix(1, nrow = n)
  if (is.null(S)) S <- matrix(1, nrow = n)

  k1 <- ncol(X); k2 <- ncol(Z); k3 <- ncol(S)

  link.mu <- stats::make.link(link$link.mu)
  link.sigma <- stats::make.link(link$link.sigma)

  if (is.null(m_start)){
    m_start <- m_inits(data, link, flambda, gen)
  }

  #### Copula ------------------------------------------------------------------
  if (is.null(P_start)){
    P_start <- cormat$start(y)
  }

  if (is.null(c_start)){
    c_start <- rep(4, ell(copula)$nextraP)
  }

  init <- c(m_start, P_start, c_start)

  ## Log-likelihood ------------------------------------------------------------
  ll <- function(par){

    ### Marginal setting -------------------------------------------------------
    beta <- as.vector(par[1:k1])
    gamma <- as.vector(par[1:k2 + k1])
    zeta <- NULL
    nu <- NULL

    if (!flambda){
      zeta <-  as.vector(par[1:k3 + k1 + k2])
      nMpar <- k1 + k2 + k3
    } else{
      zeta <- as.vector(0)
      nMpar <- k1 + k2
    }

    mu <- as.numeric(link.mu$linkinv(X%*%beta))
    sigma <- as.numeric(link.sigma$linkinv(Z%*%gamma))
    lambda <- as.numeric(S%*%zeta)
    if (BCSgen(gen)$extraP){
      nu <- par[nMpar + 1]
      nMpar <- nMpar + 1
    }

    ### Copula setting ---------------------------------------------------------
    nPpar <- cormat$npar
    nCpar <- ell(copula)$nextraP

    alpha <- NULL
    delta <- NULL

    if (nPpar > 0)  alpha <- par[1:nPpar + nMpar]
    if (nCpar > 0)  delta <- par[1:nCpar + (nMpar + nPpar)]

    P <- cormat$P(alpha, n)

    if (any(!is.finite(mu)) | any(!is.finite(sigma)) | any(!is.finite(lambda)) |
        any(mu < 0) | any(sigma < 0) | any(alpha < -1) | any(alpha > 1)){
      NaN
    }else {
      copulaGEN <- ell(copula)$gen
      epsilon <- BCSgen(copulaGEN)$q(pBCS(y, mu, sigma, lambda, nu, gen),
                                     nu = delta)
        log(ell(copula)$Md(epsilon, P, delta = delta)) -
        sum(log(BCSgen(copulaGEN)$d(epsilon^2, nu = delta))) +
        sum(log(dBCS(y, mu, sigma, lambda, nu, gen)))

    }

  }

  ## Estimates -----------------------------------------------------------------
  opt <- suppressWarnings(stats::optim(par = init,
                                       fn = ll,
                                       method = method,
                                       control = control,
                                       hessian = hessian))

  if (opt$convergence > 0)
    warning(cat("optimization failed to converge\n optim message:\n",
                opt$message))

  opt$start <- init
  opt
}