R/simul_fun_generalized_2d.R

In STCCGEV: Conditional Copula Model for Crop Yield Forecasting

#' A Special Case of simulation_generalized in 2 Dimensions
#'
#' @param nsim Integer, number of simulation replications.
#' @param n_train Integer, number of training observations.
#' @param n_test Integer, number of test observations.
#' @param copula Character, specifying the copula type: "Clayton", "Frank",
#' "Gumbel", "Joe", or "Gaussian".
#' @param init_params Numeric vector, initial parameter values for optimization.
#' @param fn Function, log-likelihood function for parameter estimation.
#' @param u1 Numeric vector (n_train), first pseudo-observation for the copula.
#' @param u2 Numeric vector (n_train), second pseudo-observation for the copula.
#' @param z1_train Numeric matrix (n_train x M), observed data for the first margin.
#' @param z2_train Numeric matrix (n_train x M), observed data for the second margin.
#' @param X_t Numeric matrix (n_train x M), risk factors for the dynamic copula parameter.
#' @param y1_test Numeric vector (n_test), true future values for the first response variable.
#' @param y2_test Numeric vector (n_test), true future values for the second response variable.
#' @param BSTS_1 Fitted BSTS model for the first response variable.
#' @param BSTS_2 Fitted BSTS model for the second response variable.
#'
#' @return A list containing:
#' \item{theta_simulated}{Simulated copula parameters across replications.}
#' \item{y1_simulated}{Simulated values for the first response variable.}
#' \item{y2_simulated}{Simulated values for the second response variable.}
#' \item{MSE}{Mean squared error for each simulation run.}
#' \item{optim_results}{Results from the optimization process.}
#'
#' @importFrom stats cor optim rnorm dnorm qnorm predict
#' @importFrom evd rgev dgev
#' @importFrom copula claytonCopula frankCopula gumbelCopula joeCopula normalCopula rCopula
#' @import bsts
#' @export
#'
## Function to optimize the full pseudo-loglikelihood and perform new forecasts
simul_fun_generalized_2d <- function(nsim,
                                     n_train,
                                     n_test,
                                     copula,
                                     init_params,
                                     fn,
                                     u1,
                                     u2,
                                     z1_train,
                                     z2_train,
                                     X_t,
                                     y1_test,
                                     y2_test,
                                     BSTS_1,
                                     BSTS_2){

  results <- list()

  if(copula != "Clayton" & copula != "Frank" & copula != "Gumbel" & copula != "Joe" & copula != "Gaussian"){stop("Invalid copula provided")}

  optim_results <- optim(par = init_params,
                         fn = fn,
                         u1 = u1,
                         u2 = u2,
                         X_t = X_t,
                         z1 = z1_train,
                         z2 = z2_train,
                         copula = copula)

  omega <- optim_results$par[1]
  alpha <- optim_results$par[2]
  gamma <- optim_results$par[3:(3+ncol(z1_train)-1)]

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

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

  #if(copula == "Gaussian"){
  #  theta_simulated <- lapply(1:6, function(x) {
  #    matrix(NA, nrow = nsim, ncol = n_test)
  #  })
  #} else{
  #theta_simulated <- matrix(NA, nrow = nsim, ncol = n_test)
  #}

  theta_simulated <- matrix(NA, nrow = nsim, ncol = n_test)

  y1_simulated <- matrix(NA, nrow = nsim, ncol = n_test)
  y2_simulated <- matrix(NA, nrow = nsim, ncol = n_test)

  z1_simulated <- matrix(NA, nrow = nsim, ncol = n_test)
  z2_simulated <- matrix(NA, nrow = nsim, ncol = n_test)

  MSE <- rep(NA, nsim)

  mu1_sim <- vector("list", ncol(z1_train))
  for(i in seq_along(mu1_sim)){
    mu1_sim[[i]] <- matrix(0, nrow = nsim, ncol = n_test)
  }

  mu2_sim <- vector("list", ncol(z2_train))
  for(i in seq_along(mu2_sim)){
    mu2_sim[[i]] <- matrix(0, nrow = nsim, ncol = n_test)
  }

  z1_sim <- vector("list", ncol(z1_train))
  for(i in seq_along(z1_sim)){
    z1_sim[[i]] <- matrix(0, nrow = nsim, ncol = n_test)
  }

  z2_sim <- vector("list", ncol(z2_train))
  for(i in seq_along(z2_sim)){
    z2_sim[[i]] <- matrix(0, nrow = nsim, ncol = n_test)
  }

  for(j in 1:ncol(z1_train)){
    mu1_sim[[j]][,1] <- colMeans(z1_train)[j]
    for(i in 1:nsim){
      z1_sim[[j]][i,1] <- rgev(1, loc = mu1_sim[[j]][i,1], scale = sigma_gev1[j], shape = xi_gev1[j])
    }
  }

  for(j in 1:ncol(z2_train)){
    mu2_sim[[j]][,1] <- colMeans(z2_train)[j]
    for(i in 1:nsim){
      z2_sim[[j]][i,1] <- rgev(1, loc = mu2_sim[[j]][i,1], scale = sigma_gev2[j], shape = xi_gev2[j])
    }
  }

  if(copula == "Clayton"){
    theta_simulated[,1] <- clayton.theta(mean(cor(cbind(u1, u2), method = "kendall")))
  }
  else if(copula == "Frank"){
    theta_simulated[,1] <- frank.theta(mean(cor(cbind(u1, u2), method = "kendall")))
  }
  else if(copula == "Gumbel"){
    theta_simulated[,1] <- GH.theta(mean(cor(cbind(u1, u2), method = "kendall")))
  }
  else if(copula == "Joe"){
    theta_simulated[,1] <- joe.theta(mean(cor(cbind(u1, u2), method = "kendall")))
  }
  else if(copula == "Gaussian"){
    theta_simulated[,1]  <- cor(u1, u2, method = "pearson")
  }

  for(i in 1:nsim){

    for(j in 1:ncol(z1_train)){

      for(t in 2:n_test){

        mu1_sim[[j]][i,t] <- phi_gev1[j] * mu1_sim[[j]][i,t-1] + rnorm(1, 0, sigma_mu1[j])
        z1_sim[[j]][i,t] <- rgev(1, loc = mu1_sim[[j]][i,t], scale = sigma_gev1[j], shape = xi_gev1[j])

        mu2_sim[[j]][i,t] <- phi_gev2[j] * mu2_sim[[j]][i,t-1] + rnorm(1, 0, sigma_mu2[j])
        z2_sim[[j]][i,t] <- rgev(1, loc = mu2_sim[[j]][i,t], scale = sigma_gev2[j], shape = xi_gev2[j])

        if(copula == "Clayton"){
          theta_simulated[i,t] <- abs(omega + alpha*theta_simulated[i,t-1] + sum(gamma*max(z1_sim[[j]][i,t-1],z2_sim[[j]][i,t-1]) ) )
        }
        else if(copula == "Frank"){
          theta_simulated[i,t] <- abs(omega + alpha*theta_simulated[i,t-1] + sum(gamma*max(z1_sim[[j]][i,t-1],z2_sim[[j]][i,t-1]) ) )
        }
        else if(copula == "Gumbel"){
          theta_simulated[i,t] <- abs(omega + alpha*theta_simulated[i,t-1] + sum(gamma*max(z1_sim[[j]][i,t-1],z2_sim[[j]][i,t-1]) ) ) + 1
        }
        else if(copula == "Joe"){
          theta_simulated[i,t] <- abs(omega + alpha*theta_simulated[i,t-1] + sum(gamma*max(z1_sim[[j]][i,t-1],z2_sim[[j]][i,t-1]) ) ) + 1
        }
        else if(copula == "Gaussian"){
          theta_simulated[i,t] <- tanh(omega + alpha*theta_simulated[i,t-1] + sum(gamma*max(z1_sim[[j]][i,t-1],z2_sim[[j]][i,t-1]) ) )
        }

      }
    }

    #u_sim <- matrix(NA, nrow = n_test, ncol = 2)
    u_sim <- matrix(NA, nrow = BSTS_1$niter, ncol = 2*n_test)

    col_index_u1 <- seq(1, 2*n_test-1, by = 2)
    col_index_u2 <- seq(2, 2*n_test, by = 2)


    #for (j in 1:nrow(u_sim)){
    for (j in 1:n_test){

      if(copula == "Clayton"){
        u_sim[,c(col_index_u1[j],col_index_u2[j])] <- rCopula(BSTS_1$niter, claytonCopula(param = theta_simulated[i,j], dim = 2))
      }
      else if(copula == "Frank"){
        u_sim[,c(col_index_u1[j],col_index_u2[j])] <- rCopula(BSTS_1$niter, frankCopula(param = theta_simulated[i,j], dim = 2))
      }
      else if(copula == "Gumbel"){
        u_sim[,c(col_index_u1[j],col_index_u2[j])] <- rCopula(BSTS_1$niter, gumbelCopula(param = theta_simulated[i,j], dim = 2))
      }
      else if(copula == "Joe"){
        u_sim[,c(col_index_u1[j],col_index_u2[j])] <- rCopula(BSTS_1$niter, joeCopula(param = theta_simulated[i,j], dim = 2))
      }
      else if(copula == "Gaussian"){
        u_sim[,c(col_index_u1[j],col_index_u2[j])] <- rCopula(BSTS_1$niter, normalCopula(param = theta_simulated[i,j], dim = 2, dispstr = "un"))
      }

    }

    # Extract the new simulates values of the covariates
    z1_new <- t(do.call(rbind, lapply(z1_sim, function(mat) mat[i, ])))
    z2_new <- t(do.call(rbind, lapply(z2_sim, function(mat) mat[i, ])))

    epsilon_y1 <- matrix(NA, nrow = BSTS_1$niter, ncol = n_test)
    epsilon_y2 <- matrix(NA, nrow = BSTS_2$niter, ncol = n_test)

    selected_columns_1 <- seq(1, 2*n_test-1, by = 2)
    selected_columns_2 <- seq(2, 2*n_test, by = 2)

    for(j in 1:n_test){
      epsilon_y1[,j] <- qnorm(u_sim[,selected_columns_1[j]], mean = 0, sd = BSTS_1$sigma.obs)
      epsilon_y2[,j] <- qnorm(u_sim[,selected_columns_2[j]], mean = 0, sd = BSTS_2$sigma.obs)
    }

    predictions_y1 <- matrix(NA, nrow = BSTS_1$niter, ncol = n_test)
    predictions_y2 <- matrix(NA, nrow = BSTS_2$niter, ncol = n_test)

    predy1 <- predict(BSTS_1, horizon = 1, newdata = z1_new, burn = BSTS_1$niter*0.1)$distribution
    predy2 <- predict(BSTS_2, horizon = 1, newdata = z2_new, burn = BSTS_2$niter*0.1)$distribution

    predictions_y1 <- apply(predy1, 2, function(col, sds) {
      col - rnorm(1, mean = 0, sd = sds)
    }, sds = BSTS_1$sigma.obs)

    predictions_y2 <- apply(predy2, 2, function(col, sds) {
      col - rnorm(1, mean = 0, sd = sds)
    }, sds = BSTS_2$sigma.obs)

    y1_simulated[i,] <- (predictions_y1 + epsilon_y1[c(BSTS_1$niter*0.1+1:BSTS_1$niter),])[sample(1),]
    y2_simulated[i,] <- (predictions_y2 + epsilon_y2[c(BSTS_2$niter*0.1+1:BSTS_2$niter),])[sample(1),]

    MSE[i] <- sum((y1_simulated[i,] - y1_test)^2 + (y2_simulated[i,] - y2_test)^2) /(2*n_test)

    # print(i/nsim)

  }


  results[[1]] <- theta_simulated
  results[[2]] <- y1_simulated
  results[[3]] <- y2_simulated
  results[[4]] <- MSE
  results[[5]] <- optim_results


  return(results)

}