R/log_likelihood_generalized_2d.R

In STCCGEV: Conditional Copula Model for Crop Yield Forecasting

#' Generalized Log-Likelihood Function for 2D Copula-GEV Model
#'
#' Computes the negative log-likelihood of a 2-dimensional copula-GEV model,
#' incorporating dynamic Generalized Extreme Value (GEV) parameters and a
#' time-varying copula structure.
#'
#' @param params Numeric vector, model parameters including copula and GEV parameters.
#' @param u1 Numeric vector (length `n_train`), pseudo-observations for margin 1.
#' @param u2 Numeric vector (length `n_train`), pseudo-observations for margin 2.
#' @param X_t Numeric matrix (`n_train x M`), risk factors affecting copula parameters.
#' @param z1 Numeric matrix (`n_train x M`), observed data for margin 1.
#' @param z2 Numeric matrix (`n_train x M`), observed data for margin 2.
#' @param copula Character, specifying the copula type: "Clayton", "Frank", "Gumbel", "Joe", or "Gaussian".
#' @return The negative log-likelihood value for optimization.
#' @importFrom stats cor rnorm dnorm
#' @importFrom evd dgev
#' @importFrom copula claytonCopula frankCopula gumbelCopula joeCopula normalCopula dCopula
#' @export
#'
#' @examples test_ll_2d <-log_likelihood_generalized_2d(init_params_full,
#'                                   uu[,1],
#'                                   uu[,2],
#'                                   xx_train,
#'                                   zz_train[,1,],
#'                                   zz_train[,2,],
#'                                   "Gaussian")
#'
log_likelihood_generalized_2d <- function(params, u1, u2, X_t, z1, z2, copula) {

  result <- try({# Parameters
    omega <- params[1]
    alpha <- params[2]
    gamma <- params[3:(3+ncol(z1)-1)]

    phi_gev1 <- params[(3+ncol(z1)):(3+(ncol(z1)*2)-1)]  # AR(1) coefficient for GEV 1
    sigma_mu1 <- params[(3+(ncol(z1)*2)):(3+(ncol(z1)*3)-1)]  # Std dev of innovations for AR(1) process for GEV 1
    sigma_gev1 <- params[(3+(ncol(z1)*3)):(3+(ncol(z1)*4)-1)]  # GEV scale parameter for GEV 1
    xi_gev1 <- params[(3+(ncol(z1)*4)):(3+(ncol(z1)*5)-1)] # GEV shape parameter for GEV 1

    phi_gev2 <- params[(3+(ncol(z1)*5)):(3+(ncol(z1)*6)-1)]   # AR(1) coefficient for GEV 2
    sigma_mu2 <- params[(3+(ncol(z1)*6)):(3+(ncol(z1)*7)-1)]  # Std dev of innovations for AR(1) process for GEV 2
    sigma_gev2 <- params[(3+(ncol(z1)*7)):(3+(ncol(z1)*8)-1)] # GEV scale parameter for GEV 2
    xi_gev2 <- params[(3+(ncol(z1)*8)):(3+(ncol(z1)*9)-1)]  # GEV shape parameter for GEV 2

    n <- length(u1)
    theta <- numeric(n)
    #if(copula == "Gaussian"){
    #  theta <- matrix(NA, nrow = n, ncol = 6)
    #} else{
    #  theta <- numeric(n)
    #}
    mu_gev1 <- matrix(NA, nrow = n, ncol = ncol(z1))  # Dynamic GEV location parameter for z1
    mu_gev2 <-  matrix(NA, nrow = n, ncol = ncol(z2))  # Dynamic GEV location parameter for z2


    # Initialize the dynamic GEV location parameters (starting values)
    mu_gev1[1,] <- colMeans(z1)
    mu_gev2[1,] <- colMeans(z2)

    X_12 <- mapply(pmax, as.data.frame(z1), as.data.frame(z2))

    # Initial theta based on copula type
    if (copula == "Clayton") {
      #theta[1] <- abs(omega + alpha * clayton.theta(cor(u1, u2, method = "kendall")) +  sum(as.numeric(gamma) * colMeans(X_t)))
      theta[1] <- abs(omega + alpha * clayton.theta(mean(cor(cbind(u1, u2), method = "kendall"))) +  sum(as.numeric(gamma) * colMeans(X_t)))
    } else if (copula == "Frank") {
      #theta[1] <- abs(omega + alpha * frank.theta(mean(cor(cbind(u1, u2), method = "kendall"))) + gamma * mean(X_t))
      theta[1] <- abs(omega + alpha * 3.58 +  sum(as.numeric(gamma) * colMeans(X_t)))
    } else if (copula == "Gumbel") {
      theta[1] <- abs(omega + alpha * GH.theta(mean(cor(cbind(u1, u2), method = "kendall"))) +  sum(as.numeric(gamma) * colMeans(X_t))) + 1
    } else if (copula == "Joe") {
      #theta[1] <- abs(omega + alpha * joe.theta(mean(cor(cbind(u1, u2), method = "kendall"))) + gamma * mean(X_t)) + 1
      theta[1] <- abs(omega + alpha * 2.003 +  sum(as.numeric(gamma) * colMeans(X_t)))+ 1
    } else if (copula == "Gaussian") {
      theta[1] <- tanh(omega + alpha * cor(u1, u2, method = "pearson") +  sum(as.numeric(gamma) * colMeans(X_12)))
    } else {
      stop("Copula provided is not valid")
    }

    # Initialize log-likelihood
    ll <- 0

    # Loop over time
    for (t in 2:n) {
      # Update dynamic GEV location parameters using AR(1)

      epsilon_mu1 <- rep(NA, ncol(z1))
      ar1_density_1 <- rep(NA, ncol(z1))
      for (j in 1:length(epsilon_mu1)){
        epsilon_mu1[j] <- rnorm(1, mean = 0, sd = sigma_mu1[j])
        mu_gev1[t,j] <- phi_gev1[j] * mu_gev1[t - 1,j] + epsilon_mu1[j]
        ar1_density_1[j] <- dnorm(epsilon_mu1[j], mean = 0, sd = sigma_mu1[j], log = TRUE)
      }
      epsilon_mu2 <- rep(NA, ncol(z2))
      ar1_density_2 <- rep(NA, ncol(z2))
      for (j in 1:length(epsilon_mu2)){
        epsilon_mu2[j] <- rnorm(1, mean = 0, sd = sigma_mu2[j])
        mu_gev2[t,j] <- phi_gev2[j] * mu_gev2[t - 1,j] + epsilon_mu2[j]
        ar1_density_2[j] <- dnorm(epsilon_mu2[j], mean = 0, sd = sigma_mu2[j], log = TRUE)
      }

      # Update conditional copula parameter theta
      if (copula == "Clayton") {
        theta[t] <- dynamic.theta.clayton(params = c(omega, alpha, gamma), theta[t - 1], X_t[t,])
        cop <- claytonCopula(param = theta[t], dim = 2)
      } else if (copula == "Frank") {
        theta[t] <- dynamic.theta.frank(params = c(omega, alpha, gamma), theta[t - 1], X_t[t,])
        cop <- frankCopula(param = theta[t], dim = 2)
      } else if (copula == "Gumbel") {
        theta[t] <- dynamic.theta.gumbel(params = c(omega, alpha, gamma), theta[t - 1], X_t[t,])
        cop <- gumbelCopula(param = theta[t], dim = 2)
      } else if (copula == "Joe") {
        theta[t] <- dynamic.theta.joe(params = c(omega, alpha, gamma), theta[t - 1], X_t[t,])
        cop <- joeCopula(param = theta[t], dim = 2)
      } else if (copula == "Gaussian") {
        theta[t] <- dynamic.rho(params = c(omega, alpha, gamma), theta[t - 1], X_12[t,])
        cop <- normalCopula(param = theta[t], dim = 2, dispstr = "un")
      }

      # GEV density contributions for z1 and z2
      gev_density_1 <- rep(NA, ncol(z1))
      gev_density_2 <- rep(NA, ncol(z2))
      for(j in 1:length(gev_density_1)){
        gev_density_1[j] <- dgev(z1[t,j], loc = mu_gev1[t,j], scale = sigma_gev1[j], shape = xi_gev1[j], log = TRUE)
        gev_density_2[j] <- dgev(z2[t,j], loc = mu_gev2[t,j], scale = sigma_gev2[j], shape = xi_gev2[j], log = TRUE)
      }

      # Copula contribution
      copula_density <- log(dCopula(cbind(u1[t], u2[t]), copula = cop))

      # Update log-likelihood
      ll <- ll + copula_density + sum(gev_density_1) + sum(gev_density_2) + sum(ar1_density_1) + sum(ar1_density_2)
    }

    # Print params and log-likelihood for debugging
    # print(c(params, ll))

    # Return negative log-likelihood for optimization
    return(-ll)
  },silent = TRUE)

  # Handle errors
  if (inherits(result, "try-error")) {
    message("Error encountered in log-likelihood calculation: ", result)
    return(Inf)  # Return a penalty value to continue optimization
  }

  return(result)
}