R/log_likelihood_Generalized.R

In STCCGEV: Conditional Copula Model for Crop Yield Forecasting

#' Compute Log-Likelihood for a Generalized Dynamic Copula-GEV Model
#'
#' Computes the log-likelihood for a time-varying copula model combined with
#' Generalized Extreme Value (GEV) margins.
#'
#' @param params Numeric vector of model parameters, including copula parameters
#' (omega, alpha, gamma) and GEV distribution parameters.
#' @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 evd dgev
#' @importFrom copula claytonCopula frankCopula gumbelCopula joeCopula normalCopula dCopula
#' @importFrom stats cor rnorm dnorm
#' @importFrom utils combn
#' @export
#'
#' @examples test_ll <- log_likelihood_Generalized(init_params_full_G,uu,
#'                       zz_train,xx_train,"Gaussian")
#'
log_likelihood_Generalized <- 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)]
    offset <- 2 + M   # everything after this belongs to GEV blocks

    #####################################################################
    # 3) Parse GEV parameters from 'params'
    #     For each (d,m), we have 4 parameters:
    #       1) phi_gev[d,m],
    #       2) sigma_mu[d,m],
    #       3) sigma_gev[d,m],
    #       4) xi_gev[d,m].
    # total = 4 * (D*M)
    #####################################################################
    phi_gev   <- matrix(NA, nrow = D, ncol = M)
    sigma_mu  <- matrix(NA, nrow = D, ncol = M)
    sigma_gev <- matrix(NA, nrow = D, ncol = M)
    xi_gev    <- matrix(NA, nrow = D, ncol = M)

    for (d_idx in seq_len(D)) {
      for (m_idx in seq_len(M)) {
        phi_gev[d_idx,   m_idx] <- params[offset + 1]
        sigma_mu[d_idx,  m_idx] <- params[offset + 2]
        sigma_gev[d_idx, m_idx] <- params[offset + 3]
        xi_gev[d_idx,    m_idx] <- params[offset + 4]
        offset <- offset + 4
      }
    }

    #####################################################################
    # 4) Initialize dynamic GEV location parameters mu_gev[t,d,m]
    #####################################################################
    mu_gev <- array(NA, dim = c(T, D, M))
    # For t=1, set each mu_gev[1,d,m] to some initial value
    for (d_idx in seq_len(D)) {
      mu_gev[1, d_idx, ] <- colMeans(Z[, d_idx,])
    }

    # We'll store the AR(1) innovation eps and densities as well
    epsilon_mu <- array(NA, dim = c(T, D, M))
    ar1_density <- array(NA, dim = c(T, D, M))  # log-density

    # We'll also store the GEV log-density for each (t,d,m)
    gev_density <- array(NA, dim = c(T, D, M))

    # For t=1, compute GEV density at the initial location
    for (d_idx in seq_len(D)) {
      for (m_idx in seq_len(M)) {
        gev_density[1, d_idx, m_idx] <- evd::dgev(
          x     = Z[1, d_idx, m_idx],
          loc   = mu_gev[1, d_idx, m_idx],
          scale = sigma_gev[d_idx, m_idx],
          shape = xi_gev[d_idx, m_idx],
          log   = TRUE
        )
      }
    }

    #####################################################################
    # 5) 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)
    }

    #####################################################################
    # 6) 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

    #####################################################################
    # 7) The main time loop: t=2..T
    #####################################################################
    for (t_idx in 2:T) {
      # (a) AR(1) update for GEV location in each (d_idx, m_idx)
      for (d_idx in seq_len(D)) {
        for (m_idx in seq_len(M)) {
          # update mu_gev[t_idx,d_idx,m_idx]
          eps_val <- stats::rnorm(1, mean=0, sd=max(sigma_mu[d_idx,m_idx],1e-6))
          mu_gev[t_idx,d_idx,m_idx] <- phi_gev[d_idx,m_idx]*mu_gev[t_idx-1,d_idx,m_idx] + eps_val

          epsilon_mu[t_idx,d_idx,m_idx] <- eps_val
          ar1_density[t_idx,d_idx,m_idx] <- stats::dnorm(
            eps_val,
            mean=0,
            sd=max(sigma_mu[d_idx,m_idx],0),
            log=TRUE
          )

          # GEV log-density
          val_gev <- evd::dgev(
            x = Z[t_idx,d_idx,m_idx],
            loc = mu_gev[t_idx,d_idx,m_idx],
            scale = sigma_gev[d_idx,m_idx],
            shape = xi_gev[d_idx,m_idx],
            log=TRUE
          )
          # if (!is.finite(val_gev)) {
          #   cat("GEV DENSITY IS -Inf or NaN at time=", t_idx,
          #       " margin=", d_idx, " sub=", m_idx, "\n")
          #   cat("   data Z=", Z[t_idx, d_idx, m_idx],
          #       "   loc=", mu_gev[t_idx, d_idx, m_idx],
          #       "   scale=", sigma_gev[d_idx, m_idx],
          #       "   shape=", xi_gev[d_idx, m_idx], "\n")
          #   # Possibly: browser() to inspect further
          # }
          gev_density[t_idx,d_idx,m_idx] <- val_gev
        }
      }

      # (b) Build the copula object & evaluate log-density
      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))

      # cat("Time t=", t_idx, " CopulaDensity=", copula_density,
      #     " sumGEV=", sum(gev_density[t_idx, , ]), "\n")
      # if (!is.finite(copula_density) || !is.finite(sum(gev_density[t_idx, , ]))) {
      #   cat("NON-FINITE DETECTED at time:", t_idx, "\n")
      #   browser()
      # }

      # (d) sum the densities for time t_idx
      ll_t <- copula_density +
        sum(ar1_density[t_idx, , ]) +
        sum(gev_density[t_idx, , ])
      ll <- ll + ll_t
    }

    # Return negative log-likelihood
    return(-ll)

  }, silent=TRUE)

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