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

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

## -----------------------------------------------------------------------------
sim <- simulate_var_example(T = 200, seed = 42)

Y <- sim$Y
head(Y)

## -----------------------------------------------------------------------------
impact_true <- sim$impact_true
dimnames(impact_true) <- list(sim$var_names, sim$shock_names)
round(impact_true, 3)

## -----------------------------------------------------------------------------
p   <- 1L
dat <- build_yX(Y, p = p)

length(dat$y_list)
sapply(dat$X_list, dim)

## -----------------------------------------------------------------------------
s2out <- build_s2(dat$y_list)
round(s2out$s2, 3)

## -----------------------------------------------------------------------------
opt <- optim_hyper(
    y_list  = dat$y_list,
    X_list  = dat$X_list,
    p       = p,
    s2      = s2out$s2
)

opt$kappa1
opt$kappa2
opt$log_ml
opt$convergence

## -----------------------------------------------------------------------------
set.seed(20260514)
mcmc <- post_sim(
    Y     = Y,
    p     = p,
    R     = 1000,
    seed  = 20260514
)

names(mcmc)
mcmc$spec$kappa
mcmc$spec$log_ml

## -----------------------------------------------------------------------------
S <- rbind(
    supply   = c( 1, -1,  0),
    demand   = c( 1,  1,  0),
    monetary = c(-1, -1,  1)
)
colnames(S) <- sim$var_names
S

## -----------------------------------------------------------------------------
accepted <- sign_restrict(mcmc, S, n_id = 3L, max_tries = 1000L, verbose = FALSE)
length(accepted)

## -----------------------------------------------------------------------------
n_tries <- attr(accepted, "n_tries")
summary(n_tries)

## -----------------------------------------------------------------------------
S_partial <- S[c("supply", "demand"), ]
S_partial

accepted_partial <- sign_restrict(mcmc, S_partial, max_tries = 1000L, verbose = FALSE)
length(accepted_partial)

## -----------------------------------------------------------------------------
accepted_h3 <- sign_restrict(mcmc, S, horizons = 3L, max_tries = 5000L, verbose = FALSE)
length(accepted_h3)
summary(attr(accepted_h3, "n_tries"))

## -----------------------------------------------------------------------------
S_array <- array(0, dim = c(3, 3, 3),
                 dimnames = list(sim$shock_names, sim$var_names, paste0("h", 0:2)))
S_array[, , "h0"] <- S
S_array["monetary", "interest_rate", "h1"] <- 1
S_array["supply", "output", "h2"] <- 1
S_array

accepted_arr <- sign_restrict(mcmc, S_array, horizons = 2L, max_tries = 5000L, verbose = FALSE)
length(accepted_arr)

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

## -----------------------------------------------------------------------------
mcmc$spec$convergence

## -----------------------------------------------------------------------------
attr(accepted, "failed_draws")

## -----------------------------------------------------------------------------
max(attr(accepted, "n_tries"))

## -----------------------------------------------------------------------------
H    <- 12L
irfs <- collect_irfs(accepted, n = 3L, p = p, H = H)
dim(irfs)

## ----fig.width = 7.5, fig.height = 6.5----------------------------------------
plot_irfs(
    irfs,
    var_names   = sim$var_names,
    shock_names = sim$shock_names,
    probs       = c(0.16, 0.84)
)

## -----------------------------------------------------------------------------
h <- 4L
resp <- irfs[, h + 1, , ]                  # variable x shock x draw
med  <- apply(resp, c(1, 2), stats::median)
lo   <- apply(resp, c(1, 2), stats::quantile, probs = 0.16)
hi   <- apply(resp, c(1, 2), stats::quantile, probs = 0.84)
dimnames(med) <- dimnames(lo) <- dimnames(hi) <- list(sim$var_names, sim$shock_names)

round(med, 3)                               # posterior median at horizon 4
round(lo, 3)                                # 16th percentile
round(hi, 3)                                # 84th percentile

## -----------------------------------------------------------------------------
path <- apply(irfs[1, , 3, ], 1, stats::quantile, probs = c(0.16, 0.5, 0.84))
colnames(path) <- paste0("h", 0:H)
round(path, 3)                              # output after a monetary shock

## -----------------------------------------------------------------------------
impact_post <- Reduce("+", lapply(accepted, `[[`, "impact")) / length(accepted)
dimnames(impact_post) <- list(sim$var_names, sim$shock_names)
round(impact_post, 3)

## -----------------------------------------------------------------------------
round(impact_true, 3)

## ----error = TRUE-------------------------------------------------------------
try({
post_sim(as.data.frame(Y), p = 1L, R = 100L)
})

## ----error = TRUE-------------------------------------------------------------
try({
post_sim(Y, p = 1.5, R = 100L)
})

## ----error = TRUE-------------------------------------------------------------
try({
sign_restrict(mcmc, S[, 1:2])
})

## ----error = TRUE-------------------------------------------------------------
try({
S_bad <- S; S_bad[1, 1] <- 2
sign_restrict(mcmc, S_bad)
})

## ----error = TRUE-------------------------------------------------------------
try({
collect_irfs(accepted, n = 2L, p = p, H = 12L)
})