R/plot_forecast_compare.R

In STCCGEV: Conditional Copula Model for Crop Yield Forecasting

#' Compare Forecasts from Two Models
#'
#' Generates a time series plot comparing the forecasts from two models along with observed data.
#'
#' @param forecast1 Numeric matrix, forecasted values from the first model (columns: time points).
#' @param forecast2 Numeric matrix, forecasted values from the second model (columns: time points).
#' @param data_train Numeric vector, training data used for modeling.
#' @param data_test Numeric vector, actual test data for evaluation.
#' @param time Numeric vector, representing the time points corresponding to the data.
#' @param quant_high Numeric, upper quantile (e.g., 0.9) for confidence interval.
#' @param quant_low Numeric, lower quantile (e.g., 0.1) for confidence interval.
#' @param col1 Character, color for observed data lines.
#' @param title Character, title for the plot.
#' @return A `ggplot2` object showing the forecast comparison.
#' @importFrom ggplot2 ggplot aes geom_line geom_ribbon scale_color_manual theme_light labs
#' @importFrom grDevices adjustcolor
#' @importFrom stats quantile
#' @export
plot_forecast_compare <- function(forecast1, forecast2, data_train, data_test, time,
                                  quant_high, quant_low, col1, title){

  # Prepare data for the first model
  Observed <- c(data_train, rep(NA, length(data_test)))
  Future <- c(rep(NA, (length(data_train)-1)), data_train[length(data_train)], data_test)
  AvgForecast1 <- c(rep(NA, (length(data_train)-1)), data_train[length(data_train)], colMeans(forecast1))
  QuantHigh1 <- c(rep(NA, length(data_train)), apply(forecast1, 2, quantile, probs=quant_high))
  QuantLow1 <- c(rep(NA, length(data_train)), apply(forecast1, 2, quantile, probs=quant_low))

  # Prepare data for the second model
  AvgForecast2 <- c(rep(NA, (length(data_train)-1)), data_train[length(data_train)], colMeans(forecast2))
  QuantHigh2 <- c(rep(NA, length(data_train)), apply(forecast2, 2, quantile, probs=quant_high))
  QuantLow2 <- c(rep(NA, length(data_train)), apply(forecast2, 2, quantile, probs=quant_low))

  Time <- as.numeric(time)
  res_data <- data.frame(Time, Observed, Future, AvgForecast1, AvgForecast2,
                         QuantHigh1, QuantLow1, QuantHigh2, QuantLow2)

  # Create the ggplot
  ts_plot <- ggplot2::ggplot(res_data, aes(x = Time)) +
    # Plot observed data
    ggplot2::geom_line(aes(y = Observed, color = "Observed data"), linewidth = 0.8, na.rm = TRUE) +
    # Plot actual forecasted data
    ggplot2::geom_line(aes(y = Future, color = "Observed data"), linetype = "dotted", linewidth = 1,na.rm = TRUE) +

    # Plot average forecast for model 1
    ggplot2::geom_line(aes(y = AvgForecast1, color = "Copula-based forecast"), size = 1,na.rm = TRUE) +
    # Plot confidence interval for model 1
    ggplot2::geom_ribbon(aes(ymin = QuantLow1, ymax = QuantHigh1), fill = adjustcolor("#CF9FFF", alpha.f = 0.25), alpha = 0.35,
                         na.rm = TRUE) +

    # Plot average forecast for model 2
    ggplot2::geom_line(aes(y = AvgForecast2, color = "Inla-based forecast"), size = 1,na.rm = TRUE) +
    # Plot confidence interval for model 2
    ggplot2::geom_ribbon(aes(ymin = QuantLow2, ymax = QuantHigh2), fill = adjustcolor("darkgrey", alpha.f = 0.25), alpha = 0.35,
                         na.rm = TRUE) +

    # Labels and theme adjustments
    ggplot2::labs(title = title,
         x = "Time", y = "Crop yields") +

    # Define colors for the lines
    ggplot2::scale_color_manual(name = "Legend",
                       values = c("Observed data" = col1,
                                  "Observed data" = col1,
                                  "Copula-based forecast" = "#CF9FFF",
                                  "Inla-based forecast" = "darkgrey")) +

    ggplot2::theme_light() +

    # Position the legend
    ggplot2::theme(legend.position = "right")

  return(ts_plot)
}