r - Replicate BGVAR plot with ggplot2

Question

I would like to replicate this plot generated in the package BGVAR with ggplot2.

Here is some information on the plot command of BGVAR: https://github.com/mboeck11/BGVAR/blob/master/R/plot.R (more specifically, check plot.bgvar.irf).

From the vignette, consider this example:

set.seed(1)
library(BGVAR)
data(eerData) 
model.ssvs.1<-bgvar(Data=eerData, 
                    W=W.trade0012,
                    draws=100,
                    burnin=100,
                    plag=1,
                    prior="SSVS",
                    hyperpara=NULL, 
                    SV=TRUE,
                    thin=1,
                    trend=TRUE,
                    h=0,
                    save.country.store=FALSE,
                    eigen=1.05
)

shocks<-list();shocks$var="stir";shocks$cN<-"US";shocks$ident="chol";shocks$scal=-100
irf.chol.us.mp<-irf(model.ssvs.1,shock=shocks,n.ahead=24,save.store=TRUE)
plot(irf.chol.us.mp,resp="US.y")

This code generates this plot, which I would like to replicate in ggplot2:

The data can be obtained by:

df <- plot(irf.chol.us.mp,resp="US.y")

How to achieve the same plot with ggplot2, using df?

Any help is very welcome. Thank you in advance!

EDIT

So far, I got here:

library(tidyverse)
library(ggplot2)

US.y <- plot(irf.chol.us.mp,resp="US.y") %>% 
  data.frame() %>%
  mutate(steps = seq(from = 0, to = 24))

ggplot(US.y, aes(x = steps)) +
  geom_line(aes(y = median)) +
  geom_polygon(aes(y = low25)) + 
  geom_polygon(aes(y = high75)) +
  geom_hline(yintercept = 0, linetype='dotted', col = 'red') +
  scale_x_continuous(breaks = c(0, 4, 8, 12, 16, 20, 24)) +
  theme(axis.title = element_blank())

I am struggling with the polygons.

question from:https://stackoverflow.com/questions/66054686/replicate-bgvar-plot-with-ggplot2

Answer 1

There are a couple of issues. The first problem is that the plot method from BGVAR only exports back the 25% and 75% confidence limits. Wheras the plot also has a ribbon showing the 16% and 84% confidence limits.

To obtain these additional data points, I wrote a slightly altered version of the package's plot function, that returns both these limits in a list. There might be a simpler way to obtain these using the package functions, but I'm not familiar with BGVAR.

f = function (x, ..., resp = NULL, shock.nr = 1, cumulative = FALSE) 
{
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  if (length(shock.nr) != 1) {
    stop("Please select only one shock.")
  }
  posterior <- x$posterior
  varNames <- dimnames(posterior)[[1]]
  varAll <- varNames
  cN <- unique(sapply(strsplit(varNames, ".", fixed = TRUE), 
                      function(x) x[1]))
  vars <- unique(sapply(strsplit(varNames, ".", fixed = TRUE), 
                        function(x) x[2]))
  if (!is.null(resp)) {
    resp.p <- strsplit(resp, ".", fixed = TRUE)
    resp.c <- sapply(resp.p, function(x) x[1])
    resp.v <- sapply(resp.p, function(x) x[2])
    if (!all(unique(resp.c) %in% cN)) {
      stop("Please provide country names corresponding to the ones of the 'bgvar.irf' object.")
    }
    cN <- cN[cN %in% resp.c]
    varNames <- lapply(cN, function(x) varNames[grepl(x, 
                                                      varNames)])
    if (all(!is.na(resp.v))) {
      if (!all(unlist(lapply(resp, function(r) r %in% varAll)))) {
        stop("Please provide correct variable names corresponding to the ones in the 'bgvar' object.")
      }
      varNames <- lapply(varNames, function(l) l[l %in% 
                                                   resp])
    }
    max.vars <- unlist(lapply(varNames, length))
  }
  else {
    varNames <- lapply(cN, function(cc) varAll[grepl(cc, 
                                                     varAll)])
  }
  for (cc in 1:length(cN)) {
    rows <- max.vars[cc]/2
    if (rows < 1) 
      cols <- 1
    else cols <- 2
    if (rows%%1 != 0) 
      rows <- ceiling(rows)
    if (rows%%1 != 0) 
      rows <- ceiling(rows)
    par(mar = BGVAR:::bgvar.env$mar, mfrow = c(rows, cols))
    for (kk in 1:max.vars[cc]) {
      idx <- grep(cN[cc], varAll)
      idx <- idx[varAll[idx] %in% varNames[[cc]]][kk]
      x <- posterior[idx, , shock.nr, c("low25", "median", 
                                        "high75"), drop = TRUE]
      y <- posterior[idx, , shock.nr, c("low16", "median", 
                                        "high84"), drop = TRUE]
      if (cumulative) {
        x <- apply(x, 2, cumsum)
        y <- apply(y, 2, cumsum)
      }
      b <- range(x, y)
      b1 <- b[1]
      b2 <- rev(b)[1]
      matplot(x, type = "l", col = c("black", BGVAR:::bgvar.env$plot$col.50, 
                                     "black"), yaxt = "n", xaxt = "n", lwd = BGVAR:::bgvar.env$plot$lwd.line, 
              ylab = "", main = varAll[idx], cex.main = BGVAR:::bgvar.env$plot$cex.main, 
              cex.axis = BGVAR:::bgvar.env$plot$cex.axis, cex.lab = BGVAR:::bgvar.env$plot$cex.lab, 
              lty = c(0, 1, 0), ylim = c(b1, b2))
      polygon(c(1:nrow(y), rev(1:nrow(y))), c(y[, 1], rev(y[, 
                                                            3])), col = BGVAR:::bgvar.env$plot$col.75, border = NA)
      polygon(c(1:nrow(x), rev(1:nrow(x))), c(x[, 1], rev(x[, 
                                                            3])), col = BGVAR:::bgvar.env$plot$col.68, border = NA)
      segments(x0 = 1, y0 = 0, x1 = nrow(x), y1 = 0, col = BGVAR:::bgvar.env$plot$col.zero, 
               lty = BGVAR:::bgvar.env$plot$lty.zero, lwd = BGVAR:::bgvar.env$plot$lwd.zero)
      lines(x[, 2], col = BGVAR:::bgvar.env$plot$col.50, lwd = 4)
      axis(2, at = seq(b1, b2, length.out = 5), labels = format(seq(b1, 
                                                                    b2, length.out = 5), digits = 1, nsmall = 1), 
           cex.axis = 1.2, las = 1)
      axisindex <- seq(1, nrow(x), by = 4)
      axis(side = 1, las = 1, at = axisindex, labels = c(0:nrow(x))[axisindex], 
           cex.axis = 1.6, tick = FALSE)
      abline(v = axisindex, col = BGVAR:::bgvar.env$plot$col.tick, 
             lty = BGVAR:::bgvar.env$plot$lty.tick)
    }
  }
  return(invisible(list(x,y)))
}

Now that we have these in a list, we can convert to a data.frame and plot using geom_ribbon (rather than geom_polygon which you had)

df1 <- f(irf.chol.us.mp,resp="US.y")

df1 <- do.call(cbind.data.frame, df1)[,-5]
df1$steps = seq_len(NROW(df1))

ggplot(df1, aes(x = steps)) +
  geom_ribbon(aes(ymin=low16, ymax=high84), fill='grey80') +
  geom_ribbon(aes(ymin=low25, ymax=high75), fill='grey50') +
  geom_line(aes(y = median)) +
  geom_hline(yintercept = 0, linetype='dotted', col = 'red') +
  scale_x_continuous(breaks = c(0, 4, 8, 12, 16, 20, 24)) +
  theme(axis.title = element_blank())