R/plot.sim.R

In bplloyd/Core:

plot.sim = function(sim, actual, quantiles=NULL, levels=NULL, ylim = NULL, lwd=2,...){
  if(!("ts" %in% class(sim))){
    t0_sim = zoo::as.Date(zoo::as.yearmon(row.names(sim)[1]))
    t0_act = zoo::as.Date(zoo::as.yearmon(row.names(actual)[1]))
    sim = ts(data = sim, start = c(lubridate::year(t0_sim), lubridate::month(t0_sim)), frequency = 12)
    actual = ts(data = actual, start = c(lubridate::year(t0_act), lubridate::month(t0_act)), frequency = 12)
  }

  #get ts time data
  tsp_act = tsp(actual)
  tsp_sim = tsp(sim)
  #calculate quantile probs from levels
  if(is.null(quantiles)){
    probs.upper = 0.5 + levels/200
    probs.lower = 0.5 - levels/200
    probs = c(probs.lower, probs.upper)
    probs = probs[order(probs)]
    quantiles=ts(t(apply(sim, 1, function(r)quantile(r, probs = probs))), start = tsp_sim[1], frequency = tsp_sim[3])
    #sim_mean = rowMeans(sim)
  } else {
    levels = 100 - 2*sapply(strsplit(colnames(quantiles)[1:2], "%"), function(s)as.numeric(s[1]))
    levels = levels[order(levels)]
  }
  #determine length of simulation and number of sims
  if(is.null(dim(sim))){
    nr = length(sim)
    nc = 1
  } else {
    nr = nrow(sim)
    nc = ncol(sim)
  }

  #create ts of actual data with times for sims appended as NAs
  ts_actual_alldates = ts(data = c(as.vector(actual), rep(NA, nr)), start = tsp_act[1], frequency = tsp_act[3])
  #actual_alldates = c(as.vector(actual), rep(NA, nr))



  #get last observation and last time in order to fill the gap
  lastObs = actual[length(actual)]
  lastTime = time(actual)[length(actual)]


  #calculate quantiles and mean of the simulated data
  #sim_quantiles = ts(t(apply(sim, 1, function(r)quantile(r, probs = probs))), start = tsp_sim[1], frequency = tsp_sim[3])

  #sim_mean = ts(rowMeans(sim), start = tsp_sim[1], frequency = tsp_sim[3])


  #append the last observation to the beginning
  # sim_quantiles = rbind(rep(lastObs, ncol(quantiles)), quantiles)
  sim_quantiles = ts(data = rbind(rep(lastObs, ncol(quantiles)), quantiles),
                     start = lastTime,
                     frequency = tsp_sim[3])
  # sim_mean = c(lastObs, rowMeans(sim))
  sim_mean = ts(data = c(lastObs,rowMeans(sim)),
                start = lastTime,
                frequency = tsp_sim[3])



  #colors for prediction intervals
  # Using very small confidence levels.
  if(min(levels) < 50){
    shadecols <- rev(colorspace::sequential_hcl(100)[levels])
  }
  # This should happen almost all the time. Colors mapped to levels.
  else {
    shadecols <- rev(colorspace::sequential_hcl(52)[levels-49])
  }

  xxx <- time(sim_quantiles)
  nint <- ncol(sim_quantiles)/2

  if(is.null(ylim)){
    rng = range(sim_quantiles)
    ylim = c(0, round(rng[2], 3) + 0.01)
  }

  x_t = time(ts_actual_alldates)
  at = which((x_t %% 1) == 0)



  #plot actual
 # plot(ts_actual_alldates, ylim = ylim, lwd=1, xaxt = "n")
  plot(ts_actual_alldates, ylim = ylim, lwd=1, xaxt = "n", ...)
  #axis(1, at = at, labels = round(x_t[at], 0))
  axis(1, at =  round(x_t[at], 0), labels = round(x_t[at], 0))


  #plot PIs
  for(i in 1:nint) {
    vec = c(sim_quantiles[, i], rev(sim_quantiles[, 2*nint - i + 1]))
    polygon(c(xxx,rev(xxx)), vec,
            col=shadecols[i], border=FALSE)
  }
  #plot mean
  lines(sim_mean, col = "blue", lwd = lwd)
}

# levels = seq(10, 90, 10)
#
#
#
#
# sim_hyb = dist_sim$Simulations$hybrid
# sim_ar = dist_sim$Simulations$individual$arima
# sim_stlm  = dist_sim$Simulations$individual$stlm
# #sim_nets = dist_sim$Simulations$individual$nnetar
#
# act = dist_sim$Actual
#
#
# if(is.null(dim(sim))){
#   nr = length(sim)
#   nc = 1
# } else {
#   nr = nrow(sim)
#   nc = ncol(sim)
# }
#
# #ts_actual_alldates = ts(data = rbind(replicate(expr = as.vector(act), n = nc), replicate(expr = rep(NA, nr), n=nc)), start = tsp_act[1], frequency = tsp_act[3])
# ts_actual_alldates = ts(data = c(as.vector(act), rep(NA, nr)), start = tsp_act[1], frequency = tsp_act[3])
#
# lastObs = act[length(act)]
# lastTime = time(act)[length(act)]
#
# sim_quantiles = ts(t(apply(sim, 1, function(r)quantile(r, probs = probs))), start = tsp_sim[1], frequency = tsp_sim[3])
# sim_mean = ts(rowMeans(sim), start = tsp_sim[1], frequency = tsp_sim[3])
#
# sim_quantiles = ts(data = rbind(rep(lastObs, ncol(sim_quantiles)), sim_quantiles),
#                    start = lastTime,
#                    frequency = tsp_sim[3])
#
# sim_mean = ts(data = c(lastObs,sim_mean),
#                    start = lastTime,
#                    frequency = tsp_sim[3])
#
#
# # Using very small confidence levels.
# if(min(x$level) < 50){
#   shadecols <- rev(colorspace::sequential_hcl(100)[levels])
# }
# # This should happen almost all the time. Colors mapped to levels.
# else {
#   shadecols <- rev(colorspace::sequential_hcl(52)[levels-49])
# }
#
# #xxx <- time(x$upper)
# xxx <- time(sim_quantiles)
#
# #idx <- rev(order(x$level))
# #idx <- rev(order(x$level))
#
# #nint <- length(x$level)
# nint <- length(levels)
#
#
# plot(ts_actual_alldates, ylim = c(0, .04))
# for(i in 1:nint) {
#   vec = c(sim_quantiles[, i], rev(sim_quantiles[, 2*nint - i + 1]))
#   polygon(c(xxx,rev(xxx)), vec,
#           col=shadecols[i], border=FALSE)
# }
#
# # else if(shaded)
# # {
# #
# #   # polygon(c(xxx,rev(xxx)), c(x$lower[,idx[i]],rev(x$upper[,idx[i]])),
# #   #         col=shadecols[i], border=FALSE)
# # }
#
#
#
#
#
#
#
#
#
# # for(i in 1:ncol(sim)){
# #   lines(sim[, i], col = "blue")
# # }