## ----setup, include = FALSE---------------------------------------------------
has_bsvarSIGNs <- requireNamespace("bsvarSIGNs", quietly = TRUE)
knitr::opts_chunk$set(
    collapse = TRUE,
    comment  = "#>",
    fig.width  = 7,
    fig.height = 5,
    dpi = 100,
    eval = has_bsvarSIGNs
)
options(mcmc.pb = FALSE)

## -----------------------------------------------------------------------------
library(acpbvar)

## -----------------------------------------------------------------------------
data("monetary", package = "bsvarSIGNs")
Y <- as.matrix(monetary)
Y[, 1:5] <- 100 * Y[, 1:5]
var_names <- c("real GDP", "GDP deflator", "commodity prices",
               "total reserves", "non-borrowed reserves", "fed funds")
colnames(Y) <- var_names
dim(Y)
round(head(Y, 3), 2)

## -----------------------------------------------------------------------------
p <- 12L
t_post <- system.time(
    mcmc <- post_sim(Y, p = p, R = 1000L, var_names = var_names, seed = 1)
)
round(mcmc$spec$kappa, 4)
mcmc$spec$convergence
t_post[["elapsed"]]

## -----------------------------------------------------------------------------
S <- matrix(c(0, -1, -1, 0, -1, 1), nrow = 1,
            dimnames = list("monetary", var_names))
S
t_sign <- system.time(
    accepted <- sign_restrict(mcmc, S, horizons = 5L, max_tries = 20000L, verbose = FALSE)
)
length(accepted)
attr(accepted, "failed_draws")
summary(attr(accepted, "n_tries"))
t_sign[["elapsed"]]

## ----fig.height = 9-----------------------------------------------------------
H <- 24L
irfs <- collect_irfs(accepted, n = 6L, p = p, H = H)
irf_mp <- irfs[, , 1, , drop = FALSE]          # the identified shock only
plot_irfs(irf_mp, var_names = var_names, shock_names = "monetary policy",
          probs = c(0.16, 0.84))

## -----------------------------------------------------------------------------
path <- apply(irf_mp[1, , 1, ], 1, stats::quantile, probs = c(0.16, 0.5, 0.84))
colnames(path) <- paste0("h", 0:H)
round(path[, c(1, 7, 13, 25)], 3)             # real GDP after a monetary tightening

## -----------------------------------------------------------------------------
sign_irf <- matrix(NA, 6, 6)
sign_irf[, 1] <- c(NA, -1, -1, NA, -1, 1)      # variable x shock, horizon 0
sign_irf <- array(sign_irf, dim = c(6, 6, 6))  # horizons 0..5
set.seed(1)
t_bs <- system.time({
    spec_bs <- bsvarSIGNs::specify_bsvarSIGN$new(Y, p = p, sign_irf = sign_irf, mc.cores = 1)
    post_bs <- bsvars::estimate(spec_bs, S = 1000, show_progress = FALSE)
    irf_bs  <- bsvars::compute_impulse_responses(post_bs, horizon = H)
})
dim(irf_bs)                                   # variable x shock x horizon x draw
t_bs[["elapsed"]]

## ----fig.height = 5.5---------------------------------------------------------
med_ac <- apply(irf_mp[, , 1, ], c(1, 2), stats::median)
med_bs <- apply(irf_bs[, 1, , ], c(1, 2), stats::median)
op <- par(mfrow = c(2, 3), mar = c(2.5, 2.5, 2, 0.5), mgp = c(1.4, 0.4, 0))
for (j in 1:6) {
    ylim <- range(med_ac[j, ], med_bs[j, ])
    plot(0:H, med_ac[j, ], type = "l", lwd = 2, ylim = ylim, xlab = "", ylab = "",
         main = var_names[j], cex.main = 0.9)
    lines(0:H, med_bs[j, ], lwd = 2, lty = 2, col = "firebrick")
    abline(h = 0, lty = 3, col = "grey50")
    if (j == 1) legend("topright", c("acpbvar", "bsvarSIGNs"), lwd = 2,
                       lty = c(1, 2), col = c("black", "firebrick"), bty = "n", cex = 0.8)
}
par(op)

## -----------------------------------------------------------------------------
data.frame(package = c("acpbvar", "bsvarSIGNs"),
           seconds = round(c(t_post[["elapsed"]] + t_sign[["elapsed"]], t_bs[["elapsed"]]), 1))

## -----------------------------------------------------------------------------
accepted_joint <- sign_restrict(mcmc, S, horizons = 5L, method = "joint", verbose = FALSE)
length(accepted_joint)