Joint Bayesian inference about impulse responses in VAR models

Applied work typically reports a structural impulse response function as a sequence of separate, horizon-by-horizon summaries — a posterior median or mean at each horizon, bracketed by pointwise credible intervals built from the marginal posterior distribution at that horizon alone — and the paper argues this practice is not generally valid as either a point estimate or a measure of uncertainty about the impulse response vector. Point estimates and error bands should instead be derived from formal Bayesian decision theory applied to the joint posterior distribution of the full response vector theta, because questions applied researchers actually care about — whether a response is hump-shaped, whether there is delayed overshooting, whether a shock generates comovements consistent with stagflation — are properties of the whole vector and its cross-horizon, cross-variable dependence, not of any single coordinate (Section 1, pp. 457-459; Section 2.2.2, pp. 461-462).

The Bayes estimator theta-hat, defined as the vector minimizing posterior expected loss, is by construction always contained in the feasible response space Theta = g(Lambda); the vector of pointwise posterior medians or means is not, in general, contained in Theta whenever the number of estimated impulse responses n_irf exceeds the number of structural parameters n_p — the typical applied case. In that case Theta is a lower-dimensional smooth manifold inside R^{n_irf} rather than the full Euclidean space, so treating each horizon’s marginal posterior median or mean as if it were unconstrained can produce an impulse-response shape that does not correspond to any possible parameterization of the underlying autoregressive model (Section 2.2, pp. 460-462).

For y_t = rho*y_{t-1} + u_t with posterior rho ~ N(0.7, 1/5), the vector of pointwise posterior means at horizons 1 through 4 is [0.70, 0.69, 0.76, 0.95] — a response that initially falls and then rises — which is a shape no single value of rho in (0,1) can generate, since the true impulse response rho^h is monotonic in h for any fixed rho. The paper’s Figure 1 shows the same pathology for the pointwise posterior median relative to the (smooth, feasible) Bayes estimator under absolute loss (Section 2.2.1-2.2.2, pp. 460-462).

The Bayes estimator solves theta-hat = argmin over theta-bar in Theta of the posterior expectation of a loss function L(theta, theta-bar), and the paper derives it under four loss functions: quadratic loss (sum of squared coordinate errors), absolute loss (sum of absolute coordinate errors), Dirac delta loss (which reduces expected-loss minimization to maximizing the joint posterior density, i.e., finding the posterior mode), and an angular loss based on the average cosine-distance between response vectors, which is invariant to the scaling of the data and of the structural shocks. Under quadratic and absolute loss the estimator is approximated by Monte Carlo as the draw-weighted minimizer over the set of M posterior draws (Eqs. 1-5; Section 2.1-2.2, Section 4.1, pp. 460-462, 470).

*The 100(1-alpha)% joint credible set is defined as the smallest-posterior-risk region {theta-bar in Theta : E[L(theta, theta-bar)] <= c_(1-alpha)} containing 100(1-alpha)% posterior probability, approximated by sorting the M posterior draws by their average loss and retaining the lowest-loss (1-alpha)100% of them — each retained draw is a full impulse-response trajectory across all horizons and variables, so plotting all of them together produces a “shot-gun trajectory plot” that preserves the shape of each response path and the comovement across variables, unlike a pointwise band built from separate horizon-by-horizon quantiles. The paper argues shot-gun plots allow direct probability statements about qualitative IRF features (hump shape, sign, delayed overshooting) that pointwise bands cannot support, since pointwise bands discard the across-horizon and across-variable dependence structure of theta (Section 2.4-2.5, pp. 463-465).

*Under Dirac delta loss the Bayes estimator is the mode of the joint posterior density of theta — the single most likely structural impulse-response vector — and the corresponding 100(1-alpha)% joint credible set is the highest-posterior-density (HPD) set, obtained by evaluating the joint posterior density at each of the M draws (via closed-form expressions the paper derives for sign-identified models, Proposition 1, Eq. 22, and for exactly identified short- and/or long-run exclusion-restriction models, Proposition 2, Eq. 31, each involving a Jacobian correction for the change of variables from reduced-form to structural parameters) and retaining the (1-alpha)100% of draws with the highest density values. These propositions generalize the authors’ earlier Dirac-delta-loss results (Inoue and Kilian 2013, 2019), which were restricted to the special case n_irf = n_p, to allow an arbitrary number of estimated impulse responses (Section 3, pp. 465-469).

In the quarterly demand/supply/monetary-policy VAR(4) for U.S. output growth, the federal funds rate, and inflation (1954:IV-2007:IV), the paper reports that “in some cases, the range spanned by the responses in the joint credible set is three times as wide” as the conventional 68% pointwise bands, illustrating that pointwise bands can substantially understate the true joint estimation uncertainty about the response vector. This distortion is present regardless of whether n_irf exceeds n_p, since it stems from ignoring cross-horizon/cross-variable dependence rather than from the manifold-feasibility problem specifically (Section 2.5, p. 464; Section 5.1, p. 472).

Across both empirical examples the absolute, quadratic, and Dirac delta loss functions yield Bayes estimates that are “quite similar” to one another, and the angular-loss estimates are likewise close to the others despite being more computationally demanding to compute, leading the authors to conclude there is “little to choose” between the loss functions for typical applications and to recommend absolute or quadratic loss as cheaper defaults. Under Dirac delta and angular loss specifically, the Bayes estimator can fail to be unique even when n_irf <= n_p; the paper’s practical fix is to report the neighborhood of draws within kappa_M = log(M)/sqrt(M) percent of the minimum expected loss, and in the empirical examples this set “consists of a singleton” in most cases and, when not, “typically only two solutions… that are quite similar” (Section 4.1-4.2, pp. 470-471; Section 5.1, p. 472; Section 6, p. 474).

The authors flag that the angular-loss estimator is more computationally costly than absolute or quadratic loss because it requires a larger number of posterior draws M, and that Dirac delta loss under partial identification requires numerical integration, which is why they recommend absolute or quadratic loss as defaults; they also note that two-dimensional “joint error bands” are simpler to plot than shot-gun trajectories but lose the shape- and comovement-information that motivates the joint approach in the first place, so they recommend shot-gun plots over both conventional and joint error bands. The paper explicitly frames its scope as VAR impulse responses under any identification scheme (Cholesky, sign restrictions, exclusion restrictions, or partially identified/set-identified models), while noting the same joint-inference logic could in principle be adapted to forecast error variance decompositions, historical decompositions, and sequences of structural shocks, without carrying out those extensions itself (Section 4.1, p. 470; Section 2.5, p. 465; Section 6, p. 474).