R/log_likelihood_noGEV.R

In STCCGEV: Conditional Copula Model for Crop Yield Forecasting

#' Compute Log-Likelihood for a Generalized Dynamic Copula Model without GEV covariates
#'
#' Computes the log-likelihood for a time-varying copula model.
#'
#' @param params Numeric vector of model parameters, including copula parameters
#' (omega, alpha, gamma).
#' @param U Numeric matrix (n_train x D), pseudo-observations for the copula.
#' @param Z Numeric array (n_train x D x M), observed data for each margin and sub-feature.
#' @param X Numeric matrix (n_train x M), risk factors for the dynamic copula parameter.
#' @param copula Character, specifying the copula type: "Clayton", "Frank",
#' "Gumbel", "Joe", or "Gaussian".
#' @return Numeric, negative log-likelihood value.
#' @importFrom copula claytonCopula frankCopula gumbelCopula joeCopula normalCopula dCopula
#' @importFrom stats cor rnorm dnorm
#' @importFrom utils combn
#' @export
#'
#' @examples test_ll_noGEV <- log_likelihood_noGEV(init_params_noGEV,uu,
#'                       zz_train,x_train,"Gaussian")
#'
log_likelihood_noGEV <- function(
    params,
    U,         # n_train x D matrix of pseudo-observations for the copula
    Z,         # n_train x D x M array of data for each margin & sub-feature
    X,         # n_train x M matrix of risk factors for the dynamic copula parameter
    copula     # string: "Clayton", "Frank", "Gumbel", "Joe", "Gaussian"
    ){
      # Safely wrap in try(...) to handle numeric issues
      out <- try({

        #####################################################################
        # 1) Dimensions
        #####################################################################
        T <- dim(Z)[1]   # number of time points
        D <- dim(Z)[2]   # number of margins (for the copula)
        M <- dim(Z)[3]   # number of sub-features for each margin

        # Check U dimension
        if (nrow(U) != T || ncol(U) != D) {
          stop("U must be (T x D) to match T and number of copula margins D.")
        }

        #####################################################################
        # 2) Parse copula parameters from 'params'
        #####################################################################

        omega <- params[1]
        alpha <- params[2]
        gamma <- params[3:(2 + M)]

        #####################################################################
        # 3) Copula parameter storage: theta[t]
        #####################################################################
        if (copula == "Gaussian"&&D>2) {
          # Potentially store pairwise correlations in each row
          nPairs <- D*(D-1)/2
          theta_mat <- matrix(NA, nrow = T, ncol = nPairs)
        } else {
          theta_vec <- numeric(T)
        }

        #####################################################################
        # 4) Initialize theta at t=1
        #####################################################################
        if (copula == "Clayton") {
          raw_val <- omega + alpha*clayton.theta(mean(stats::cor(U, method = "kendall"))) + sum(gamma * colMeans(X))
          theta_vec[1] <- abs(raw_val)

        } else if (copula == "Frank") {
          raw_val <- omega + alpha*3.58 + sum(gamma * colMeans(X))
          theta_vec[1] <- abs(raw_val)

        } else if (copula == "Gumbel") {
          raw_val <- omega + alpha*GH.theta(mean(stats::cor(U, method = "kendall"))) + sum(gamma * colMeans(X))
          theta_vec[1] <- abs(raw_val) + 1

        } else if (copula == "Joe") {
          raw_val <- omega + alpha*2.003 + sum(gamma * colMeans(X))
          theta_vec[1] <- abs(raw_val) + 1

        } else if (copula == "Gaussian" &&D>2) {
          pair_idx <- combn(D, 2)
          nPairs   <- ncol(pair_idx)
          for (pp in seq_len(nPairs)) {
            i <- pair_idx[1, pp]
            j <- pair_idx[2, pp]
            r_ij   <- stats::cor(U[, i], U[, j], method="pearson")
            rawval <- omega + alpha*r_ij + sum(gamma*colMeans(pmax(Z[,i,],Z[,j,])))
            theta_mat[1, pp] <- tanh(rawval)
          }
        } else if (copula == "Gaussian" && D==2) {
          raw_val <- omega+alpha*cor(U[,1],U[,2], method= "pearson") +sum(gamma*colMeans(pmax(Z[,1,],Z[,2,])))
          theta_vec[1] <-tanh(raw_val)
        } else {
          stop("Unsupported copula type.")
        }

        # We'll define an accumulator for the log-likelihood
        ll <- 0

        #####################################################################
        # 5) The main time loop: t=2..T
        #####################################################################
        for (t_idx in 2:T) {

          if (copula == "Clayton") {
            theta_vec[t_idx] <- dynamic.theta.clayton(params = c(omega, alpha, gamma), theta_vec[t_idx-1], X[t_idx,])
            cop <- copula::claytonCopula(param=theta_vec[t_idx], dim=D)
          } else if (copula == "Frank") {
            theta_vec[t_idx] <- dynamic.theta.frank(params = c(omega, alpha, gamma), theta_vec[t_idx-1], X[t_idx,])
            cop <- copula::frankCopula(param=theta_vec[t_idx], dim=D)
          } else if (copula == "Gumbel") {
            theta_vec[t_idx] <- dynamic.theta.gumbel(params = c(omega, alpha, gamma), theta_vec[t_idx-1], X[t_idx,])
            cop <- copula::gumbelCopula(param=theta_vec[t_idx], dim=D)
          } else if (copula == "Joe") {
            theta_vec[t_idx] <- dynamic.theta.joe(params = c(omega, alpha, gamma), theta_vec[t_idx-1], X[t_idx,])
            cop <- copula::joeCopula(param=theta_vec[t_idx], dim=D)
          } else if (copula == "Gaussian" &&D>2) {
            pair_idx <- combn(D, 2)
            nPairs <- ncol(pair_idx)
            R <- diag(D)
            for (pp in seq_len(nPairs)) {
              i <- pair_idx[1, pp]
              j <- pair_idx[2, pp]
              param_lag <- theta_mat[t_idx-1, pp]
              theta_mat[t_idx, pp] <- dynamic.rho(params= c(omega, alpha, gamma), theta_mat[t_idx-1,pp], pmax(Z[t_idx,i,],Z[t_idx,j,]))
              R[i,j] <- theta_mat[t_idx, pp]
              R[j,i] <- theta_mat[t_idx, pp]
            }
            cop <- copula::normalCopula(param = R[upper.tri(R)], dim=D, dispstr="un")
          } else if (copula =="Gaussian"&&D==2){
            theta_vec[t_idx] <- dynamic.rho(params = c(omega, alpha, gamma), theta_vec[t_idx-1], pmax(Z[t_idx,1,],Z[t_idx,2,]))
            cop <-copula::normalCopula(param=theta_vec[t_idx], dim =2, dispstr = "un")
          }

          # Evaluate copula density at U[t_idx, ]
          copula_density <- log(copula::dCopula(U[t_idx, ], cop))
          ll <- ll + copula_density
        }

        # Return negative log-likelihood
        return(-ll)

      }, silent=TRUE)

      # error handling
      if (inherits(out, "try-error")) {
        message("Error in copula log-likelihood: ", out)
        return(Inf)
      }
      return(out)
    }