acpbvar: a complete workflow

This vignette walks through the full estimation pipeline of acpbvar on a small simulated VAR for which the data-generating process is known. The objective is twofold: to illustrate the public interface end to end, and to verify that the package recovers a known structural impact matrix from a sign-restriction identification.

1. A structural VAR with a known sign pattern

The helper simulate_var_example() returns a simulated sample of three variables (output, inflation, interest rate) generated by a structural VAR(1) with three structural shocks (supply, demand, monetary policy). The structural impact matrix is calibrated to deliver the canonical macroeconomic sign pattern: a positive supply shock raises output and lowers inflation; a positive demand shock raises both; a positive monetary-policy shock raises the interest rate while lowering output and inflation.

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

Y <- sim$Y
head(Y)
#>         output    inflation interest_rate
#> [1,] 1.5999687  0.009553379    0.56087474
#> [2,] 1.4689501 -0.665480622    0.54289572
#> [3,] 0.7485236  0.607942746    0.08197592
#> [4,] 1.6380594 -0.570413529    0.07332373
#> [5,] 0.8958783  0.331286665   -0.40549007
#> [6,] 0.6700644  0.152115603   -0.67038980

The true structural impact matrix, against which posterior estimates will be compared, is:

impact_true <- sim$impact_true
dimnames(impact_true) <- list(sim$var_names, sim$shock_names)
round(impact_true, 3)
#>               supply demand monetary
#> output           0.5    0.2     -0.2
#> inflation       -0.6    0.2     -0.2
#> interest_rate    0.0    0.2      0.4

2. Equation-by-equation preparation

The A-form representation requires equation-by-equation response vectors and regressor matrices. The regressors include the lagged variables, an intercept, and (for equations beyond the first) the contemporaneous variables ordered before the current equation, entered with a negative sign so that the impact matrix is lower triangular by construction.

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

length(dat$y_list)
#> [1] 3
sapply(dat$X_list, dim)
#>      [,1] [,2] [,3]
#> [1,]  199  199  199
#> [2,]    4    5    6

Univariate AR fits supply the residual variances that enter the asymmetric Minnesota prior as scaling factors.

s2out <- build_s2(dat$y_list)
round(s2out$s2, 3)
#> [1] 0.347 0.379 0.184

3. Empirical Bayes hyperparameter selection

optim_hyper() selects the own-lag and cross-lag tightness parameters \((\kappa_1, \kappa_2)\) by maximising the closed-form log marginal likelihood under the asymmetric Minnesota prior. Nelder-Mead optimisation runs in log-space to guarantee positivity.

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

opt$kappa1
#> [1] 0.5417438
opt$kappa2
#> [1] 0.08853819
opt$log_ml
#> [1] -485.6492
opt$convergence
#> [1] 0

4. Posterior simulation

post_sim() wraps the full pipeline: it prepares the data, fits the AR variances, runs Empirical Bayes for the tightness, builds the asymmetric prior, and draws R independent posterior values of the structural coefficients and variances, equation by equation (exact Normal-Inverse-Gamma draws, no Markov chain).

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

names(mcmc)
#> [1] "samples" "spec"    "dat"
mcmc$spec$kappa
#>     kappa1     kappa2 
#> 0.54174385 0.08853819
mcmc$spec$log_ml
#> [1] -485.6492

5. Structural identification through sign restrictions

The sign-restriction matrix encodes the three identifying patterns. Rows index shocks (supply, demand, monetary); columns index variables (output, inflation, interest rate). A value of +1 requires a positive response on impact, -1 a negative response, and 0 leaves the response unrestricted.

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

sign_restrict() implements the Rubio-Ramirez, Waggoner and Zha (2010) rotation algorithm: for each posterior draw of the reduced form, candidate rotations are drawn from the Haar measure and screened against all 2^n_id sign flips; the first accepted rotation is stored.

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

The diagnostic vector attached to the result records, for each accepted draw, how many Q candidates were screened before acceptance.

n_tries <- attr(accepted, "n_tries")
summary(n_tries)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    1.00    8.00   19.00   29.55   41.00  192.00

Partial identification

Partial identification means that only a subset of the structural shocks is identified. The matrix S has fewer rows than the number of variables. The columns of the rotation beyond n_id are left unrestricted and carry no label; the default n_id = nrow(S) is all that is needed.

S_partial <- S[c("supply", "demand"), ]
S_partial
#>        output inflation interest_rate
#> supply      1        -1             0
#> demand      1         1             0

accepted_partial <- sign_restrict(mcmc, S_partial, max_tries = 1000L, verbose = FALSE)
length(accepted_partial)
#> [1] 1000

If a row of zeros is passed, it is rejected: it would declare an identified shock without restrictions, whose sign would then be fixed by the order in which sign flips are tried rather than by an economic hypothesis.

Restrictions beyond impact

So far the signs have been imposed on impact only (h = 0, the default horizons = 0). The argument horizons = K extends them to horizons 0, …, K, in two possible ways. The first is the same signs at every horizon: the same S is imposed on all horizons 0, …, K, as in Uhlig (2005). The responses at the various horizons are computed for each candidate rotation, so more rotations are screened per draw.

accepted_h3 <- sign_restrict(mcmc, S, horizons = 3L, max_tries = 5000L, verbose = FALSE)
#> Warning: 1 of 1000 posterior draws dropped: no rotation found within max_tries
#> = 5000. If max(attr(., 'n_tries')) is close to max_tries, increase it.
length(accepted_h3)
#> [1] 999
summary(attr(accepted_h3, "n_tries"))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>     1.0    22.0    57.0   109.6   130.0  3291.0

The second is different signs at different horizons: S is an array of size n_id x n_var x (horizons + 1), one matrix per horizon; a slice of zeros leaves that horizon unrestricted, and each shock must carry at least one restriction somewhere. In the example below there are signs for all restricted variables at impact, then only the interest rate at h = 1 and only output at h = 2.

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
#> , , h0
#> 
#>          output inflation interest_rate
#> supply        1        -1             0
#> demand        1         1             0
#> monetary     -1        -1             1
#> 
#> , , h1
#> 
#>          output inflation interest_rate
#> supply        0         0             0
#> demand        0         0             0
#> monetary      0         0             1
#> 
#> , , h2
#> 
#>          output inflation interest_rate
#> supply        1         0             0
#> demand        0         0             0
#> monetary      0         0             0

accepted_arr <- sign_restrict(mcmc, S_array, horizons = 2L, max_tries = 5000L, verbose = FALSE)
length(accepted_arr)
#> [1] 1000

Partial identification works the same way: the matrix or array S has n_id rows, and the columns beyond n_id remain free at each horizon.

Two ways of accepting a draw

The default is method = "per_draw", which searches up to max_tries rotations for each draw from the posterior and keeps the first one that satisfies the constraints. In this case the reduced-form posterior is left untouched, and Q is uniform within the set identified by each posterior draw. The alternative is method = "joint", which draws a single rotation per sample and discards the draw if that rotation does not pass. As in the replication code of Chan (2022), the reduced-form draws are implicitly reweighted by the size of their identified set. In this case far fewer draws are accepted, so a larger R is required. Chan (2022, footnote 13) describes both implementations.

accepted_joint <- sign_restrict(mcmc, S, method = "joint", verbose = FALSE)
length(accepted_joint)
#> [1] 29
attr(accepted_joint, "method")
#> [1] "joint"

Reading the diagnostics

After a run, three numbers are worth checking before looking at any impulse response: the convergence code of the hyperparameter search, the draws that found no valid rotation, and how close the search came to its limit. They are printed in this order below.

mcmc$spec$convergence
#> [1] 0

convergence is the code returned by optim() in the search for the hyperparameters \(\kappa_1\) and \(\kappa_2\). If convergence is 0, convergence has been reached; otherwise the \(\kappa\) values stored in spec are the point at which the optimiser stopped, and post_sim() has already issued a warning.

attr(accepted, "failed_draws")
#> integer(0)

failed_draws contains the indices of the draws for which no rotation satisfied the sign constraints within max_tries. If the identified set of those draws is empty, it makes sense to discard them. The problem arises only when the limit bites before a valid rotation is found, which the next number tells.

max(attr(accepted, "n_tries"))
#> [1] 192

n_tries records, for each accepted draw, the attempt at which the rotation satisfying the signs was found; max_tries is 10000 by default. Comparing max(n_tries) with max_tries allows you to determine whether the rejected draws were rejected by the model or by the limit. If max(n_tries) is close to max_tries, it is advisable to raise the limit and run again.

6. Impulse responses

collect_irfs() stacks the per-draw impulse responses into a single four-dimensional array; plot_irfs() summarises the posterior into a median path and a credible band.

H    <- 12L
irfs <- collect_irfs(accepted, n = 3L, p = p, H = H)
dim(irfs)
#> [1]    3   13    3 1000

plot_irfs(
    irfs,
    var_names   = sim$var_names,
    shock_names = sim$shock_names,
    probs       = c(0.16, 0.84)
)

The impact responses (horizon zero) carry the sign restrictions by construction: the supply column is positive on output and negative on inflation; the demand column is positive on both; the monetary column is negative on output and inflation, positive on the interest rate.

Extracting the responses

irfs is a 4-dimensional array: variable × horizon × shock × draw. The h-th horizon is located at position h + 1 because the first element represents the impact at h = 0. A summary at a fixed horizon is obtained by slicing the second index and taking quantiles over the last one, the draws: below, the posterior median and the 16-84 band at h = 4.

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
#>               supply demand monetary
#> output         0.039  0.116   -0.051
#> inflation     -0.006  0.057   -0.021
#> interest_rate -0.043  0.099    0.079
round(lo, 3)                                # 16th percentile
#>               supply demand monetary
#> output        -0.001  0.074   -0.111
#> inflation     -0.028  0.027   -0.054
#> interest_rate -0.078  0.041    0.019
round(hi, 3)                                # 84th percentile
#>               supply demand monetary
#> output         0.082  0.157    0.007
#> inflation      0.020  0.087    0.012
#> interest_rate -0.008  0.145    0.136

By fixing the variable and the shock instead of the horizon, we obtain the entire trajectory of a response.

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
#>         h0     h1     h2     h3     h4     h5     h6     h7     h8     h9
#> 16% -0.351 -0.251 -0.188 -0.143 -0.111 -0.085 -0.067 -0.052 -0.041 -0.032
#> 50% -0.185 -0.145 -0.108 -0.075 -0.051 -0.034 -0.023 -0.015 -0.009 -0.006
#> 84% -0.050 -0.054 -0.029 -0.008  0.007  0.017  0.022  0.023  0.023  0.022
#>        h10    h11    h12
#> 16% -0.025 -0.021 -0.017
#> 50% -0.004 -0.002 -0.001
#> 84%  0.020  0.018  0.016

7. Recovery of the true impact

Averaging the impact column across accepted draws gives a posterior estimate of the structural impact matrix, which can be compared against the true matrix used in the simulation.

impact_post <- Reduce("+", lapply(accepted, `[[`, "impact")) / length(accepted)
dimnames(impact_post) <- list(sim$var_names, sim$shock_names)
round(impact_post, 3)
#>               supply demand monetary
#> output         0.468  0.177   -0.199
#> inflation     -0.514  0.190   -0.204
#> interest_rate -0.038  0.159    0.309

round(impact_true, 3)
#>               supply demand monetary
#> output           0.5    0.2     -0.2
#> inflation       -0.6    0.2     -0.2
#> interest_rate    0.0    0.2      0.4

The posterior recovers the sign pattern of the true impact matrix on every entry, and the magnitudes are in the same neighbourhood, with the expected sampling and identification uncertainty. Sign restrictions identify the impact matrix only up to an orthogonal transformation within the admissible set, so exact numerical recovery is not the relevant standard; the recovery of the qualitative pattern, and of the order of magnitude of the entries, is.

8. Common mistakes

Every exported function first checks its arguments and stops, displaying a message that identifies the incorrect argument and says what is expected. The chunks below are executed with error = TRUE, so the messages are shown instead of interrupting the vignette.

The data must be a numeric matrix, not a data frame:

post_sim(as.data.frame(Y), p = 1L, R = 100L)
#> Error:
#> ! `Y` must be a numeric matrix (T x n), one column per variable.

The lag order must be a whole number:

post_sim(Y, p = 1.5, R = 100L)
#> Error:
#> ! `p` must be a single integer with 1 <= p <= 199; got 1.5.

S must have one column per variable, and entries in {-1, 0, 1}:

sign_restrict(mcmc, S[, 1:2])
#> Error:
#> ! `S` has 2 columns but the model has 3 variables.

S_bad <- S; S_bad[1, 1] <- 2
sign_restrict(mcmc, S_bad)
#> Error:
#> ! `S` must be a numeric matrix (n_id x n_var) or a 3-d array (n_id x n_var x (horizons + 1)) with entries in {-1, 0, 1}.

n must match the number of variables in the accepted draws:

collect_irfs(accepted, n = 2L, p = p, H = 12L)
#> Error:
#> ! `n` = 2 does not match the 3 x 3 impact matrices in `accepted`.