R/plot_arimas_monte.R

In HydrosystemsGroup/Weather-Generator: Daily weather generator for performing a climate stress test.

#' Statistical Distribution Plot of Monte Carlo ARIMA Simulations
#'
#' @param obs numeric vector of historical annual precipitation
#' @param wgen_mc output of monte carlo simulation of arimas generated by arimas_monte()
#' @export
plot_arimas_monte <- function(obs, wgen_mc) {
  wgen_stats <- data.frame(MEAN=colMeans(wgen_mc$x),
                           SD=apply(wgen_mc$x, 2, sd),
                           SKEW=apply(wgen_mc$x, 2, skewness))
  wgen_stats <- dplyr::mutate(wgen_stats, TRIAL=seq(1:n()))
  wgen_stats <- tidyr::gather(wgen_stats, STAT, VALUE, MEAN:SKEW)

  hist_stats <- data.frame(STAT=c('MEAN', 'SD', 'SKEW'),
                           VALUE=c(mean(obs), sd(obs), skewness(obs)))

  p <- ggplot2::ggplot(wgen_stats, aes(x=STAT, y=VALUE)) +
    ggplot2::geom_boxplot(fill='grey90') +
    ggplot2::stat_summary(fun.y=mean, geom='point', color='black', size=4) +
    ggplot2::geom_point(data=hist_stats, color='red', size=4) +
    ggplot2::facet_wrap(~STAT, ncol=3, scales='free') +
    ggplot2::labs(y='') +
    ggplot2::theme_bw() +
    ggplot2::theme(axis.ticks.x=ggplot2::element_blank(),
                   axis.title.x=ggplot2::element_blank(),
                   axis.text.x=ggplot2::element_blank())
  p
}