Scenario Analysis with Multivariate Bayesian Machine Learning Models

bibunit\section{Introduction} In this paper, we discuss how to conduct scenario analysis with dynamic multivariate models in macroeconomics when the functional form of the conditional mean is nonlinear and/or unknown. While related tools exist in linear and traditional nonlinear frameworks (e.g., regime switching or time-varying parameter models, see fischer2023general for a recent example), they rely on potentially restrictive parametric assumptions and are not directly applicable in nonparametric models. The latter feature in many recent macroeconometric applications that employ Bayesian machine learning (ML) methods huber2022inference,clark2021tail,huber2023nowcasting,hauzenberger2024gaussian,chernis2025bayesian,lima2025minnesota. We use the term scenario analysis broadly to refer to different counterfactual experiments. These include versions of conditional forecasts (CFs) and structural scenario analysis antolin2021structural, and several variants of nonlinear impulse response functions (IRFs). While some aspects and issues in the Bayesian ML context have been discussed in isolation in the aforementioned papers, a unified or comprehensive treatment is not yet available. We bridge this gap by developing a framework and estimation algorithms for scenario analysis that are widely applicable when the conditional mean is flexibly learned from the data. Conditional forecasts simulate future paths of variables under scenarios encoded as constraints on other observables, structural shocks, or both. Related techniques have been used by academics and practitioners since the 1980s in linear models doan1984forecasting, with subsequent refinements waggoner1999conditional,andersson2010density,jarocinski2010conditional,banbura2015conditional,antolin2021structural,chernis2024decision,chan2023conditional.\footnote{For a Bayesian decision analysis perspective and discussion of additional related literature, see west2024perspectives.} Breaking the assumption of linearity, however, complicates matters: there is no general representation of nonlinear time series as functions of shocks potter2000nonlinear, and unknown nonlinearities make it difficult to derive multi-step ahead predictive distributions. Obtaining closed form solutions is often impossible, but in this case one can resort to Monte Carlo methods. We show how common approaches can be adapted to nonparametric settings. Akin to banbura2015conditional, our approach casts conditional forecasting as a nonlinear state space problem. We use Particle Gibbs with Ancestor Sampling lindsten14a, which combines generality with computational efficiency, making it suitable for a comparatively wide range of nonlinear multivariate models. An IRF can be defined as the difference between conditional forecasts gallant1993nonlinear,koop1996impulse,rambachan2021common,jordalocal. Our framework thus lends itself to estimating nonlinear dynamic effects gonccalves2021impulse,gonccalves2024state,kolesar2024dynamic in the form of generalized IRFs (GIRFs), and thereby contributes to improving interpretability of ML methods. Compared with IRFs in linear settings (which are symmetric, shape invariant and history independent) the GIRFs do not feature such potentially restrictive properties. Specifically, we consider (1) unorthogonalized GIRFs (the difference between conditional and unconditional forecasts subject to certain restrictions on observables), we explore (2) structural GIRFs to shocks identified with approaches typically used in the structural VAR (SVAR) literature, and (3) restricted GIRFs, which can be used to construct policy counterfactuals and quantify the contributions of specific sets of transmission channels in the propagation of shocks.\footnote{A related approach derived from antolin2021structural is discussed in breitenlechner2022goes. Counterfactual experiments in this spirit relate to alternative policy rules and are subject to the lucas1976 critique mckay2023can. While our baseline approach is also susceptible to these concerns, there are options to improve robustness, and nonparametric modeling offers additional flexibility and a safeguard; see Section (ref) for a more detailed discussion.} The proposed approach is generally applicable in multivariate dynamic models. This means that there are many potential candidates for estimating conditional mean functions. These choices include (but are not limited to) regression trees or Gaussian process priors marcellino2024bookchapter. Such nonparametric methods are more robust to misspecification that may arise from assuming tightly parameterized likelihoods, and can be used to capturing genuinely nonlinear patterns in the data. For our empirical work, we use Bayesian Additive Regression Trees chipman2010bart as a specific nonparametric implementation. We pick this sum-of-trees model because tree-based approaches have proven particularly capable of producing accurate forecasts when used with time series for the US economy, with datasets structured similarly to the one we use in this paper medeiros2021forecasting,goulet2022machine,clark2021tail,goulet2024bag. Our approach to predictive inference is developed under the assumption of an additive multivariate Gaussian error term with a time-varying variance covariance matrix. We illustrate our framework through several empirical applications. Our dataset comprises about $25$ quarterly macroeconomic and financial variables for the US economy ranging from the mid-$1970$s to the last quarter of $2024$. In one of our explorations, we also add some international variables for the euro area (EA) and the United Kingdom (UK). The applied work assesses and illustrates the role of nonlinearities when interest centers on CFs, and we explore asymmetries in the propagation of shocks of different types, signs and magnitudes. Specifically, we provide three empirical applications. First, inspired by chan2023conditional, we use a subset of the macroeconomic assumptions underlying the annual stress test conducted by the Federal Reserve System and compute CFs using constraints on observables for different scenarios, comparing predictive densities from linear and nonlinear models. Second, reflecting the growth-at-risk literature adrian2019vulnerable, we study the unorthogonalized counterfactual implications of varying financial conditions on tail risks of output growth, inflation, and employment. Third, we identify a US-based financial shock recursively barnichon2022effects, and compute GIRFs to shocks of different signs and magnitudes that are allowed to propagate internationally. We then use a restricted GIRF approach to gauge the role of spillovers and spillbacks, inspired by the empirical work of breitenlechner2022goes. The rest of this paper is structured as follows. Section (ref) addresses the challenges and solutions to obtain predictive inference in the presence of nonlinearities of unknown form. We discuss how to impose constraints on forecasts, and how these constraints may be used to construct scenarios through GIRFs. Section (ref) presents an econometric implementation using BART. Section (ref) provides three empirical illustrations. The last section concludes. \section{Nonlinearity and Predictive distributions} Let $\bm{y}_t = (y_{1t},\hdots,y_{nt})'$ collect $n$ variables for $t = 1,\hdots, T,$ and $\bm{x}_t = (\bm{y}_{t-1}',\hdots,\bm{y}_{t-p}')'$ is a $k=np$ vector of lags. Interest centers on dynamic multivariate models of the form: \begin{equation} \bm{y}_t = \bm{F}(\bm{x}_t) + \bm{\epsilon}_t, \quad \bm{\epsilon}_t \sim \mathcal{N}(\bm{0}_n,\bm{\Sigma}_t), \end{equation} where $\bm{F}(\bm{x}_t) = (f_1(\bm{x}_t),\hdots,f_n(\bm{x}_t))'$ is an $n$-vector of conditional mean functions $f_i(\bm{x}_t):\mathbb{R}^k\rightarrow\mathbb{R}$ for $i=1,\hdots,n,$ such that $\bm{F}(\bm{x}_t):\mathbb{R}^k\rightarrow\mathbb{R}^n$. We assume iid reduced form Gaussian errors $\bm{\epsilon}_t$, with $n\times n$ time-varying covariance matrix $\bm{\Sigma}_t$. One may assume functional forms for the $f_i(\bm{x}_t)$'s or treat them as unknown and estimate them nonparametrically. The methods we propose are designed specifically for the latter case. All scenario analyses we discuss in this paper rely on computing functions of the (conditional) moments of $\bm{y}_{t+h}$ based on ((ref)). Consider the single-lag case $p = 1$, without loss of generality. We may write: \begin{equation} \bm{y}_{t+h} = \tilde{\bm{F}}_{(h)}(\bm{y}_t,\bm{\epsilon}_{t+1},\bm{\epsilon}_{t+2},\hdots,\bm{\epsilon}_{t+h}) = \bm{F}(\tilde{\bm{F}}_{(h-1)}(\bm{y}_t,\bm{\epsilon}_{t+1},\hdots,\bm{\epsilon}_{t+(h-1)})) + \bm{\epsilon}_{t+h}, \end{equation} with $\bm{y}_{t+1} = \tilde{\bm{F}}_{(1)}(\bm{y}_t,\bm{\epsilon}_{t+1}) = \bm{F}(\bm{y}_{t}) + \bm{\epsilon}_{t+1}$, and $\tilde{\bm{F}}_{(h)}(\bullet)$ denotes the $h$-step composition of $\bm{F}(\bullet)$ which is defined recursively. This expression is related to the Wold decomposition, and expresses $\bm{y}_{t+h}$ as a function of the initial condition $\bm{y}_t$ and a sequence of white noise shocks, $\{\bm{\epsilon}_{t+1},\hdots,\bm{\epsilon}_{t+h}\}$; see also gourieroux2023nonlinear and Appendix (ref).\footnote{Under some assumptions about the behavior of $\tilde{\bm{F}}_{(h)}(\bullet)$, it is noteworthy that expansions such as the Volterra series could be used as an approximation and a close equivalent to the Wold representation for nonlinear time series, see, e.g., potter2000nonlinear and jorda2005estimation.} Our framework assumes that the functional form of $\bm{F}(\bullet)$ is unknown and estimated, so it is not necessarily additively separable. Even if it were, applying a nonlinear function to a Gaussian shock does not necessarily yield a random variable that follows a well-known distribution. Thus, it is not generally possible to derive closed form higher-order moments. In order to explore the predictive distributions nevertheless, one can resort to simulation-based methods. \subsection{Predictive Simulation} Define a vector $\bm{\Xi}$ that contains all coefficients and latent variables necessary to parameterize ((ref)). At time $\tau$, the one-step ahead predictive distribution is: \begin{equation} p_\tau(\bm{y}_{\tau+1} \:\vert\:\mathcal{I}) = \int p_\tau(\bm{y}_{\tau+1} \:\vert\:\mathcal{I}, \bm{\Xi}) p(\bm{\Xi} \:\vert\:\mathcal{I}) d \bm{\Xi}, \end{equation} where $\mathcal{I}$ denotes the information set used to infer $\bm{\Xi}$, and the subscript on $p_\tau(\bullet)$ marks the forecast origin. For typical out-of-sample exercises, $\mathcal{I} = \{\bm{y}_t\}_{t = 1}^{\tau}$, and we are interested in $p_\tau(\bm{y}_{\tau+h} \:\vert\:\{\bm{y}_t\}_{t = 1}^{\tau})$ for $h = 1, 2, \hdots,$ steps ahead geweke2010comparing. In other cases we condition on $\mathcal{I}=\{\bm{y}_t\}_{t = 1}^{T}$, to compute scenarios in-sample for $\tau \in \{1,2,\hdots, T\}$ using $p_\tau(\bm{y}_{\tau+h} \:\vert\:\{\bm{y}_t\}_{t = 1}^{T})$ with parameters updated based on the full sample, but for initial conditions at time $\tau$. In either case (unless necessary for clarity, we omit the $\tau$ index for the forecast origin below), $p_\tau(\bm{y}_{\tau+1} \:\vert\:\mathcal{I})$ generally does not take a well-known form, and neither does the distribution of higher-order forecasts for $h\geq2$. However, we may obtain random samples from them via simulation. This involves exploiting the fact that even though $p(\bm{y}_{\tau+1} \:\vert\:\mathcal{I})$ is unknown, $p(\bm{y}_{\tau+1} \:\vert\:\mathcal{I}, \bm{\Xi})$ takes a known form under the model in ((ref)). The one-step ahead predictive distribution is: \begin{equation} p(\bm{y}_{\tau+1} \:\vert\:\mathcal{I}, \bm{\Xi}^{(m)}) = \mathcal{N}(\bm{F}^{(m)}(\bm{x}_{\tau+1}), \bm{\Sigma}_{\tau+1}^{(m)}), \end{equation} where $\bm{x}_{\tau+1} = (\bm{y}_{\tau}',\hdots,\bm{y}_{\tau-p+1}')'$. Here, $x^{(m)}$ refers to the $m$th draw of a random variable $x$ --- in most cases throughout this paper, $m$ indexes a draw from a posterior or predictive distribution when running a Markov chain Monte Carlo (MCMC) algorithm. In what follows we omit the $m$-indexing for the parameters, but stress that these computations are carried out in each sweep of the algorithm and thus account for the posterior uncertainty of parameters. For $h \geq 2$ we iterate forward, conditioning recursively on the draws for preceding horizons, by setting the predictors to $\bm{x}_{\tau+h}^{(m)}=(\bm{y}_{\tau+h-1}^{(m)\prime},\bm{y}_{\tau+h-2}^{(m)\prime},\hdots)'$, and obtain: \begin{equation} p(\bm{y}_{\tau+h} \:\vert\:\mathcal{I}, \bm{y}_{\tau+1:\tau+h-1}^{(m)}, \bm{\Xi}) = \mathcal{N}(\bm{F}(\bm{x}_{\tau+h}^{(m)}), \bm{\Sigma}_{\tau+h}), \end{equation} where $\bm{y}_{\tau+1:\tau+h-1}^{(m)}$ denotes the path of the variables from $\tau+1$ to $\tau+h-1$ and $\bm{y}_{\tau+1:\tau+h} = (\bm{y}_{\tau+1}',\hdots,\bm{y}_{\tau+h}')'$. This relates to the recursive composition discussed in the context of ((ref)), and exploits the fact that: \begin{equation} p(\bm{y}_{\tau+1:\tau+h} \:\vert\:\mathcal{I}) = \int p(\bm{y}_{\tau+1} \:\vert\:\mathcal{I},\bm{\Xi}) \prod_{j = 2}^h p(\bm{y}_{\tau+j} \:\vert\:\bm{y}_{\tau+1:\tau+j-1}, \mathcal{I},\bm{\Xi}) p(\bm{\Xi} \:\vert\:\mathcal{I}) d\bm{\Xi}, \end{equation} that is, the joint predictive distribution factors into the product of the conditional one-step ahead densities. Simulating the process forward, sampling from the distribution in ((ref)) across horizons $h = 1,2,\hdots,$ in each sweep of an MCMC algorithm, delivers draws from $p(\bm{y}_{\tau+1:\tau+h} \:\vert\:\mathcal{I})$ via a Monte Carlo approach. \subsection{Conditional Forecasts} As their name suggests, CFs imply an additional conditioning argument in ((ref)). In this context, we denote by $\mathcal{C}_{h}$ a set that defines desired restrictions at horizon $h = 1, 2, \hdots,H,$ where $H$ is the maximum forecast horizon of interest, $\mathcal{C}_{1:H} = \{\mathcal{C}_{1},\hdots,\mathcal{C}_{H}\}$ collects all restrictions over the full forecast path, and the unconditional forecast results when $\mathcal{C}_{h} = \emptyset$ for all $h$. The joint distribution of interest is $p(\bm{y}_{\tau+1:\tau+H} \:\vert\:\mathcal{I}, \mathcal{C}_{1:H})$, and can be obtained by marginalizing $p(\bm{y}_{\tau+1:\tau+H} \:\vert\:\mathcal{I}, \mathcal{C}_{1:H},\bm{\Xi})$ over the parameters via MCMC sampling, as discussed previously. When there are no restrictions (i.e., for unconditional forecasts), the model can be simulated forward according to the law of motion in ((ref)), so that any jointly realized random path over $(\tau+1):(\tau+H)$ is consistent with the model's dynamics. In the presence of restrictions, however, the key challenge in the nonlinear context is that we cannot work directly with the joint distribution to impose them across all horizons simultaneously, as is possible in linear settings; and a decomposition like in ((ref)) is generally not applicable. We propose to use a particle MCMC algorithm Andrieu2010, adapted from the literature on nonlinear state space models, as a solution. Specifically, our approach is based on PGAS as proposed by lindsten14a, akin to the filtering-based conditional forecast implementations of banbura2015conditional, and capable of enforcing the restrictions discussed in more detail below on an $h$-by-$h$ basis. In the spirit of andersson2010density, we define $\mathcal{C}_{h}$ as a stochastic restriction:\footnote{The restrictions can also be written as a regression, $\bm{r}_h = \bm{R}_h\bm{y}_{\tau+h} + \bm{\eta}_{h}$ with $\bm{\eta}_{h}\sim\mathcal{N}(\bm{0},\bm{\Omega}_h)$, which is equivalent in distribution. The regression illustrates the interpretation as a type of measurement equation. Exact (hard) restrictions can be approximated by setting the respective elements of $\bm{\Omega}_h$ to, e.g., $10^{-8}$, which has computational and numerical advantages over using exact zeroes.} \begin{align} \bm{R}_h\bm{y}_{\tau+h}\sim\mathcal{N}\left(\bm{r}_h,\bm{\Omega}_h\right), \end{align} where $\bm{R}_h$ is an $r_h \times n$ horizon-specific matrix selecting or encoding $r_h$ restrictions on linear combinations of endogenous variables, and $\bm{r}_h$ and $\bm{\Omega}_h$ are the $r_h\times1$ mean and $r_h\times r_h$ covariance matrix of the restrictions. \sffamilyStructural scenarios. In many cases it is desirable to impose a mix of restrictions on future observables as well as structural shocks, for structural scenario analysis in the spirit of antolin2021structural. This allows for imposing that specific structural shocks (rather than an unrestricted combination of them) are responsible for the respective scenario in terms of the observables. Econometrically, this can be achieved by forcing a subset of the structural shocks to their unconditional distribution (non-driving shocks), and leave the remaining ones (driving shocks) free. The latter thus may deviate from their unconditional distribution, thereby acting as the offsetting forces that yield the restrictions imposed on the observables. The same intuition can be used for computing dynamic responses to specific structural shocks, see also breitenlechner2022goes. As discussed above, the one-step ahead prediction of our model conditional on any preceding information up to $\tau+h-1$ can be written as $\bm{y}_{\tau+h} = \bm{\mu}_{\tau+h}^{(m)} + \bm{\epsilon}_{\tau+h}$, with $\bm{\mu}_{\tau+h}^{(m)} = \bm{F}(\bm{y}_{\tau+h-1}^{(m)}, \hdots, \bm{y}_{\tau+h-p}^{(m)}) = \bm{F}(\bm{x}_{\tau+h}^{(m)})$ and $\bm{\epsilon}_{\tau+h} \sim \mathcal{N}(\bm{0},\bm{\Sigma}_{\tau+h})$. We abstract below from time-variation in the variances to simplify notation, but note that this framework is applicable also in the presence of heteroskedastic shocks. A structural form of the model can be obtained as $\bm{H}\bm{y}_{\tau+h} = \bm{H}\bm{\mu}_{\tau+h} + \bm{u}_{\tau+h}$, with iid shocks $\bm{u}_{\tau}\sim\mathcal{N}(\bm{0}_n,\bm{I}_n)$, where $\bm{I}_n$ is an $n$-dimensional identity matrix. Then, $\bm{\epsilon}_{\tau} = \bm{H}^{-1}\bm{u}_{\tau}$, and $\bm{\Sigma} = \bm{H}^{-1}\bm{H}^{-1}{'}$.\footnote{An alternative, see also Section (ref), is to parameterize ((ref)) directly as $\bm{H}\bm{y}_t = \bm{F}(\bm{x}_t) + \bm{u}_t,~\bm{u}_t\sim\mathcal{N}(\bm{0}_n,\bm{I}_n)$, where $\bm{H}$ is nonsingular, see also arias2023macroeconomic,chan2024large. Time-varying variances of structural shocks, or more general forms of heteroskedasticity, can be addressed by using $\bm{u}_{\tau}\sim\mathcal{N}(\bm{0}_n,\bm{S}_\tau)$, where $\bm{S}_\tau = \text{diag}(s_{1\tau}^2,\hdots,s_{n\tau}^2)$, such that $\bm{\Sigma}_\tau = \bm{H}^{-1}\bm{S}_\tau\bm{H}^{-1}{'}$.} Both types of restrictions can be implemented with distributional assumptions about the restricted vector of endogenous variables, as $r_h = r_h^{(y)} + r_h^{(u)}$ restrictions on: \begin{enumerate}[label=(\roman*)] • {observables}: $\bm{R}^{(y)}_h\bm{y}_{\tau+h} \sim \mathcal{N}(\bm{r}^{(y)}_h, \bm{\Omega}^{(y)}_h)$, imposing $r_h^{(y)}$ restrictions; • {structural shocks}: $\bm{R}^{(u)}\bm{u}_{\tau+h} \sim \mathcal{N}(\bm{r}^{(u)}_h,\bm{\Omega}^{(u)}_h)$, or equivalently, $\bm{R}^{(u)}\bm{H}\bm{y}_{\tau+h} \sim \mathcal{N}(\bm{r}^{(u)}_h + \bm{R}^{(u)}\bm{H}\bm{\mu}_{\tau+h},\bm{\Omega}^{(u)}_h)$, imposing $r_h^{(u)}$ restrictions. \end{enumerate} We may write these in terms of ((ref)) by stacking both types in terms of their implications about the observables. The operator $\text{bdiag}(\bullet)$ outputs a block diagonal matrix, with the inputs arranged along the main block diagonal; we obtain $\bm{r}_h = [\bm{r}^{(y)}_h; \bm{r}^{(u)}_h + \bm{R}^{(u)}\bm{H}\bm{\mu}_{\tau+h}]$, $\bm{R}_h = [\bm{R}^{(y)}_{h}; \bm{R}^{(u)}_{h}\bm{H}]$ and $\bm{\Omega}_h = \text{bdiag}(\bm{\Omega}^{(y)}_h, \bm{\Omega}^{(u)}_h)$. Under these assumptions we obtain the joint distribution of the forecast and the restrictions at horizon $h$: \begin{equation} \begin{bmatrix} \bm{y}_{\tau+h}\\ \bm{r}_h \end{bmatrix} \sim \mathcal{N}\left( \begin{bmatrix} \bm{I}_n\\ \bm{R}_h \end{bmatrix} \bm{\mu}_{\tau+h}^{(m)}, \begin{bmatrix} \bm{\Sigma} & \bm{\Sigma}\bm{R}_h'\\ \bm{R}_h\bm{\Sigma} & \bm{R}_h\bm{\Sigma}\bm{R}_h' + \bm{\Omega}_h \end{bmatrix} \right). \end{equation} which encodes a nonlinear state combined with a linear measurement equation. This casts conditional forecasting and structural scenario analysis as a state space problem, similar to banbura2015conditional, and allows to apply the PGAS algorithm of lindsten14a with a few adjustments.\footnote{Our framework is also applicable for nowcasting problems cimadomo2022nowcasting.} \sffamilyParticle Gibbs with Ancestor Sampling. The idea of particle Gibbs algorithms is to use hypothetical realizations, a set of particles $v = 1,\hdots,V,$ combined with a scheme to score how plausible these particles are in light of the measurements/restrictions. A key ingredient of PGAS is that we condition on a fixed trajectory for one of the particles. This reference particle is chosen as the previous draw for the respective states in the encompassing MCMC sampler, ensuring that the Markov kernel preserves the correct invariant distribution Andrieu2010. The “ancestor” part of PGAS adds a resampling step for the parents of the reference trajectory, which alleviates path degeneracy by maintaining variability in the ancestral lineages, see lindsten14a for details. This improves mixing, especially with few particles, making the algorithm computationally attractive. The $h$-specific distribution of the forecast conditional on the restrictions is available in closed form west1997bayesian when assuming additive Gaussian errors and restrictions as in ((ref)): \begin{align} \bm{y}_{\tau+h} &\:\vert\:\bm{r}_h,\bm{R}_h,\bm{\Omega}_h,\bullet \sim \mathcal{N}(\bm{r}^{\ast}_h, \bm{\Omega}_h^{\ast})\\ \bm{r}^{\ast}_h &= \bm{\mu}_{\tau+h}^{(m)} + \bm{\Sigma}\bm{R}_h'(\bm{R}_h\bm{\Sigma}\bm{R}_h' + \bm{\Omega}_h)^{-1}(\bm{r}_h - \bm{R}_h \bm{\mu}_{\tau+h}^{(m)}),\nonumber\\ \bm{\Omega}_h^{\ast} &= \bm{\Sigma} - \bm{\Sigma}\bm{R}_h' (\bm{R}_h\bm{\Sigma}\bm{R}_h' + \bm{\Omega}_h)^{-1} \bm{R}_h\bm{\Sigma}.\nonumber \end{align} While ((ref)) places the respective restriction only at a single horizon, this conditional Gaussian distribution can be used as an “optimal proposal” for particles in algorithm, in the spirit of a forward-filtering update. This allows to generate particles that satisfy the restrictions locally without being exposed to excessive weight degeneracy (the situation in which only very few, if any, particles are compatible with the restrictions). Used within an encompassing MCMC algorithm, PGAS propagates these particles in such a way that the restrictions are enforced jointly over the entire forecast path, and we obtain draws from $p(\bm{y}_{\tau+1:\tau+H} \:\vert\:\mathcal{I}, \mathcal{C}_{1:H})$. Details about PGAS and explicit computations appear in Appendix (ref). We provide a brief summary of the main steps below. The initial conditions at each forecast origin $\tau$ are known, and we can generate $V-1$ candidate particles to be considered alongside the reference particle from the previous MCMC draw. We loop through: \begin{enumerate}[label=(\roman*)] • Resampling and ancestor sampling: To only retain particles with a high probability of having generated the measurements, i.e., the restrictions in ((ref)), we sample ancestors for the non-reference particles using weights from the previous horizon, see step (iii) below. For the reference particle, we sample an ancestor index from a categorical distribution with support $\{1,\hdots,V\}$, where the weights are proportional to the likelihood of particle $v = 1,\hdots,V,$ having generated the reference. Put simply, we obtain a random ancestry that is likely to have generated the reference path. • Propagation: We then propagate the “surviving” non-reference particles from the resampling step one period forward. In case there are no restrictions ($\mathcal{C}_h = \emptyset$) this is done using ((ref)). In case there are restrictions, we use ((ref)) to draw particles already conditioned on them at $h$. This avoids generating unsuitable candidate realizations and reduces the required number of particles $V$. • \textit{Weighting}: We update the weights used in step (i) for the next horizon based on ((ref)). Specifically, we compute the likelihood of the lineage of the respective particles in light of the restrictions. In case there are no restrictions, we obtain equal weights. \end{enumerate} Using the weights at the maximum forecast horizon $H$, we draw a particle index and generate a fully smoothed trajectory by tracing back through the stored ancestor indices. This trajectory becomes the new reference path in the next iteration. These steps are run in each sweep of the main MCMC algorithm. In some cases we require the expectation of the predictive distribution in each MCMC sweep. We obtain this moment via exploiting a backward recursion godsill2004monte to obtain smoothing weights, and use the full set of particle trajectories for related computations, see also Appendix (ref). We provide a comparison of the precision sampler of chan2023conditional with our PGAS approach in a linear VAR (with closed form solutions) using artificial data in Appendix (ref). \subsection{Generalized Impulse Response Functions} Various types of IRFs are widely used for both academic and policy analysis. These types of dynamic effects can generally be defined as the difference between two forecasts, see gallant1993nonlinear,koop1996impulse. This provides a natural link to our previous discussions on predictive distributions. Indeed, due to the recursive nonlinear structure of our model, we need to resort to a simulation-based version of the GIRF.\footnote{Simulation-based methods are often used to compute (generalized) impulse response functions in nonlinear parametric baumeister2013time,alessandri2019financial, or nonparametric huber2022inference,clark2025nonparametric,hauzenberger2024gaussian models.} We consider three main variants, which are nested in the expression: \begin{equation} \bm{\Delta}_{\tau} = \mathbb{E}(\bm{y}_{\tau+1:\tau+H}\:\vert\:\mathcal{C}_{1:H}^{(\texttt{s})},\bm{x}_{\tau+1},\mathcal{I}) - \mathbb{E}(\bm{y}_{\tau+1:\tau+H}\:\vert\:\mathcal{C}_{1:H}^{(\texttt{b})},\bm{x}_{\tau+1},\mathcal{I}), \end{equation} and differentiated by distinct conditioning assumptions. In line with the terminology in crump2021large, the first expectation refers to the “scenario” (\texttt{s}) forecast imposed with the restrictions $\mathcal{C}_{1:H}^{(\texttt{s})}$, and the latter is the “baseline” (\texttt{b}) forecast subject to $\mathcal{C}_{1:H}^{(\texttt{b})}$. Note that the conditioning on $\mathcal{I}$ here refers to the fact that these quantities are computed with draws for parameters conditioning on all available information; the expectation is taken at time $\tau$, as indicated by conditioning on $\bm{x}_{\tau+1}$. For later reference, define horizon specific GIRFs $\bm{\delta}_{\tau,h}$ based on the structure of $\bm{\Delta}_{\tau} = (\bm{\delta}_{\tau,1}^{\prime},\hdots,\bm{\delta}_{\tau,H}^{\prime})'$. That is, the one-step ahead forecast horizon is associated with the \textit{contemporaneous} impact of the restrictions breitenlechner2022goes. Equation ((ref)) nests our main variants:\footnote{Other related work includes bernanke1997systematic,hamilton2004comment,sims2006does,baumeister2012unconventional,baumeister2014real,adrian2025scenario.} the (1) unorthogonalized GIRF (UGIRF) which can be obtained from a fully reduced form model, but which may be extended to a structural scenario analysis when additionally identifying and restricting structural shocks; the (2) structural GIRF (SGIRF) in response to a structurally identified shock imposed exclusively with restrictions on the contemporaneous ($h=1$) shocks; and, the (3) restricted GIRF (RGIRF), using an SGIRF as a baseline but systematically manipulating transmission channels in the spirit of assessing alternate policy rules, by additionally imposing restrictions on observables in a certain way. \textsc{\textbf{Unorthogonalized GIRF}}. This version is closely related to the original GIRF of koop1996impulse. It is obtained by restricting only observables as the scenario $\mathcal{C}_{1:H}^{(\texttt{s})}$, which can be compared to the unconditional forecast, $\mathcal{C}_{1:H}^{(\texttt{b})} = \emptyset$, as a baseline. Without any further restrictions, UGIRFs simply reflect a likely combination of structural shocks that drive the change in the respective observable(s), see also crump2021large. In case one is willing to impose additional structure in the model via identifying structural shocks and $\bm{H}$, restrictions on observables can be complemented with restrictions on shocks antolin2021structural. Specifically, one may restrict a subset of the shocks to their unconditional distribution (non-driving shocks), while the unrestricted subset (driving shocks) is allowed to deviate to deliver the structural scenario in terms of the observables. In light of the assumptions about the structural form above, this implies $\bm{r}^{(u)}_h = \bm{0}$ and $\bm{\Omega}^{(u)}_h = \bm{I}$ for the \textit{non-driving} shocks. \textsc{\textbf{Structural GIRF}}. This variant relies on orthogonalized structural economic shocks. We follow the literature jordalocal, and define the SGIRF in response to the $j$th structural shock of magnitude $d$ as: \begin{equation} \bm{\Delta}_{j\tau}^{(d)} = \mathbb{E}(\bm{y}_{\tau+1:\tau+H} \:\vert\:u_{j\tau+1} = d_0 + d, \bm{x}_{\tau+1}, \mathcal{I}) - \mathbb{E}(\bm{y}_{\tau+1:\tau+H} \:\vert\:u_{j\tau+1} = d_0, \bm{x}_{\tau+1}, \mathcal{I}). \end{equation} The parameter $d_0$ is the baseline level of the shock. Since the structural errors enter linearly in our model, choices about the baseline level do not matter and we implicitly use $d_0 = 0$ in most of our subsequent discussions. Equation ((ref)) is a special case of ((ref)), which can be obtained by defining $\mathcal{C}_1^{(\texttt{s})}$ with $\bm{r}_1^{(u)} = (d_0 + d)\cdot\bm{e}_j'$, and $\mathcal{C}_1^{(\texttt{b})}$ with $\bm{r}_1^{(u)} = d_0\cdot\bm{e}_j'$; $\bm{e}_j'$ is the $j$th column of $\bm{I}_n$. For both forecasts we use $\bm{\Omega}_1^{(u)} = \text{diag}(1,\hdots,1,10^{-8},1,\hdots,1)$ with an approximately binding restriction in the $j$th position, $\bm{R}_1^{(u)} = \bm{I}_n$, and we have $\mathcal{C}_h = \emptyset$ for $h > 1$. The GIRF at time $\tau$ and horizon $h$ in response to the $j$th structural shock of size $d$ is referred to as $\bm{\delta}_{j\tau,h}^{(d)}$. In Appendix (ref) we discuss an alternative computation method which can be used without running PGAS, and Appendix (ref) again provides an illustration and comparison in a linear context using artificial data. \textsc{\textbf{Restricted GIRF}}. The final variant we consider combines the SGIRF with restrictions on observables, so that we may impose $\mathbb{E}(\bm{R}^{(y)}_h \bm{\delta}_{\tau,h}^{(d)}) = \bm{0}_{r_h^{(y)}}$ along the desired dimensions. This approach switches off specific transmission channels of structural shocks, by partially matching the observables of the scenario forecast conditional on the shock with those of the forecasted observables in the baseline predictive distribution. breitenlechner2022goes provide a discussion of this approach in a linear context, which allows for directly restricting the IRF. In our nonlinear setting, we need to account for varying initial conditions and work with the two conditional forecasts, as we cannot directly manipulate the IRF. From an implementation perspective, we may obtain a draw $\bm{y}_{\tau+1:\tau+H}^{(\texttt{b},m)}$ from the baseline distribution $p(\bm{y}_{\tau+1:\tau+H} \:\vert\:\mathcal{C}_1^{(\texttt{b})}, \bm{x}_{\tau+1}, \mathcal{I})$, where $\mathcal{C}_1^{(\texttt{b})}$ features the impact restrictions as in ((ref)) plus any desired restrictions on the future sequence of non-driving shocks (this is however not a necessary requirement, as otherwise all shocks may be allowed to deviate from their unconditional distribution; or, one can also force all of them to follow their unconditional distribution, see below). For computing the scenario forecast, we then augment the shock impact restriction in ((ref)) with a restriction on the respective observable, so that $\bm{r}_h^{(y)} = \bm{R}_h^{(y)}\bm{y}_{\tau+h}^{(\texttt{b},m)}$, to define $\mathcal{C}_{1:H}^{(\texttt{s})}$. Individual \textit{realized paths} from the scenario and baseline forecast distributions are identical along the restricted dimensions. The restricted GIRF then arises as the difference in \textit{expected values} as in ((ref)), which has an expected value of zero across MCMC iterations along the restricted dimensions, but this cannot strictly be enforced for each individual draw. This is because otherwise the reference trajectory in each sweep of the PGAS algorithm may not necessarily respect the restrictions. \textsc{\textbf{Conditional and unconditional GIRFs}}. One may obtain all three of these GIRFs conditional on each point in time $\tau$. In principle, one may thus consider “time-varying” dynamic effects of shocks for each period individually (these estimates are also referred to as “filtered” in the sense that they condition on the realized history up to a specific period, see rambachan2021common). This time variation, however, is exclusively and mechanically due to variation across initial conditions, since the functional form in the class of models we consider is potentially nonlinear but time-invariant. Unconditional (G)IRFs can be computed in various ways (e.g., by randomizing over initial conditions and averaging, in bootstrap-type approaches); see also kilian2017structural, and gonccalves2021impulse,gonccalves2024state for definitions. In case we require estimates of unconditional versions, our preferred approach is to compute the GIRF at each point in time, and then average out this source of randomness, $\overline{\bm{\Delta}} = \sum_{\tau = 1}^T \bm{\Delta}_{\tau}/T$. It is worth noting that under several assumptions, see rambachan2021common, temporally averaged GIRFs can identify a dynamic causal effect (an “average treatment effect,” in related terminology). \textsc{\textbf{Lucas critique}}. Conditional forecasts and the GIRF variants relate to work on evaluating counterfactual policy rules, and may be exposed to the lucas1976 critique. This is the case particularly when imposing highly unusual conditioning scenarios, as opposed to mere “modest policy interventions,” which may leave behavioral patterns of economic agents unchanged leeper2003modest. A recent discussion of these aspects is also provided in mckay2023can; a related issue is that any imposed counterfactual paths on observables in this spirit require specific neutralizing \textit{future} shocks, and the restrictions hold only ex post and not necessarily in expectation. Since our framework is based on the idea of a future sequence of shocks delivering the scenario as in antolin2021structural, similar concerns as in linear models apply. Two further remarks are relevant. First, by virtue of our nonparametric approach, structural shocks may locally affect relationships among variables, even if the underlying functional relationships remain constant over time. Thus, while our baseline setting is susceptible to the lucas1976 critique, it is likely more robust to related concerns than linear models. Second, breitenlechner2024fiscal recently recognized that the lucas1976-robust approach of mckay2023can can be implemented in linear frameworks of the antolin2021structural type by allowing offsetting shocks \textit{exclusively} on impact. This restriction can straightforwardly be implemented in our setup, provided the matrix $\bm{H}$ is identified. In the spirit of our discussion in Section (ref), to subsequently explore the future conditional moments of the observables stochastically, we further augmented the desired impact restriction with a restriction forcing all future structural shocks (for $h > 1$) to adhere to their unconditional distribution. \section{Model Specification and Estimation Algorithm} \subsection{Multivariate System Estimation} To estimate the model in ((ref)) we rely on a conditional representation of its $n$ equations. Let $\bm{e}_i$ of size $1 \times n$ denote the $i$th row of $\bm{I}_n$, and $\bm{E}_i$ of size $(n-1) \times n$ results from deleting the $i$th row of $\bm{I}_n$. Using $\bm{y}_{-it} = \bm{E}_i \bm{y}_t$ we may write $\bm{y}_t = \bm{e}_i'y_{it} + \bm{E}_i'\bm{y}_{-it}$. Under the assumptions of ((ref)), $p(y_{it}\:\vert\:\bm{y}_{-it},\bullet) \propto \exp\left\{-(\bm{e}_i\bm{\Sigma}_{t}^{-1}\bm{e}_i'y_{it}^2 - 2 y_{it}\bm{e}_i\bm{\Sigma}_t^{-1}(\bm{F}(\bm{x}_t) - \bm{E}_i'\bm{y}_{-it}))/2\right\}$, which is a Gaussian with variance $\varsigma_{it}^2 = (\bm{e}_i\bm{\Sigma}_{t}^{-1}\bm{e}_i')^{-1}$ and mean ${\mu}_{it} = \varsigma_{it}^2 (\bm{e}_i\bm{\Sigma}_t^{-1}(\bm{F}(\bm{x}_t) - \bm{E}_i'\bm{y}_{-it}))$. This distribution is equivalent to a common representation of the conditional multivariate Gaussian cong2017fast. The mean can alternatively be written as ${\mu}_{it} = f_i(\bm{x}_t) - \varsigma_{it}^2(\bm{e}_i\bm{\Sigma}_{t}^{-1}\bm{E}_i')(\bm{y}_{-it} - \bm{E}_i \bm{F}(\bm{x}_t))$, and we define $\tilde{\mu}_{it} = -\varsigma_{it}^2(\bm{e}_i\bm{\Sigma}_{t}^{-1}\bm{E}_i')(\bm{y}_{-it} - \bm{E}_i \bm{F}(\bm{x}_t))$, i.e., $\mu_{it} = f_i(\bm{x}_t) + \tilde{\mu}_{it}$.\footnote{Related papers often either use a mapping between the structural and reduced form of the VAR to enable equation-by-equation estimation hauzenberger2024gaussian, or rely on factor models for the reduced form errors clark2021tail. These approaches come with computational and inferential advantages and disadvantages. The former is simple to implement but requires parameterizing a structural form, which may cause issues such as inadvertently (instead of purposefully to achieve structural identification) breaking order-invariance of the equations. The latter allows for order-invariant inference but gives rise to the usual identification challenges of factor models. Our approach uses an order-invariant reduced form model for estimating the conditional mean relationships arias2023macroeconomic,chan2024large.} The $i$th equation of the multivariate model in regression form, conditional on all other equations, is then given by: \begin{equation} (y_{it} - \tilde{\mu}_{it}) = f_i(\bm{x}_t) + u_{it}, \quad u_{it} \sim \mathcal{N}(0,\varsigma_{it}^2), \end{equation} which can be used in a Gibbs sampler to update the conditional mean relationships by looping through equations $i = 1,\hdots,n$. This approach is similar to the one of esser2024seemingly and allows to treat each equation of the multivariate system individually, conditional on all other equations, which is computationally quick. \subsection{Bayesian Additive Regression Trees} The approach we discuss in Section (ref) works with any implementation of multivariate models with additive and jointly Gaussian errors. That is, assuming a linear functional form for $\bm{F}(\bm{x}_t)$ combined with suitable priors results in a standard Bayesian VAR (BVAR).\footnote{When we consider linear versions of our model for comparisons, we implement this setting with $\bm{F}(\bm{x}_t) = \bm{A}\bm{x}_t$ where $\bm{A}$ is an $n\times k$ matrix of reduced form VAR coefficients. We assume a horseshoe prior with a single global shrinkage component on these parameters, see also hauzenberger2024bookchapter; ((ref)) can be used to update the VAR coefficients equation-by-equation from textbook Gaussian posteriors.} In case we treat $\bm{F}(\bm{x}_t)$ nonparametrically, several options are available. Due to its versatility and established favorable empirical properties we mentioned earlier, we use BART to approximate the equation-specific functions in our applied work. That is, we consider a sum of $s = 1,\hdots,S,$ tree functions $\ell_{is}(\bm{x}_t\:\vert\:\mathcal{T}_{is}, \bm{\mathrm{m}}_{is})$ such that $f_i(\bm{x}_t) \approx \sum_{s=1}^{S}\ell_{is}(\bm{x}_t\:\vert\:\mathcal{T}_{is}, \bm{\mathrm{m}}_{is})$ where $\mathcal{T}_{is}$ are regression trees and $\bm{\mathrm{m}}_{is}$ is a vector of terminal node parameters (which serve as fitted values). Instead of having a single but complex tree, BART is akin to ensemble methods, and uses a sum of many simple trees (“weak learners”), which has been shown to work well. Using BART requires an algorithm that estimates splitting variables and thresholds for which we specify suitable priors that together yield $p(\mathcal{T}_{is})$; we further need a prior on the terminal node parameters $p(\bm{\mathrm{m}}_{is}\:\vert\:\mathcal{T}_{is})$. Our setup follows chipman2010bart and we first define the probability that a tree ends at a specific node at depth $d = 0,1,2,\hdots,$ as $\alpha/(1+d)^\beta$, with $\alpha\in(0,1)$ and $\beta\in\mathbb{R}^{+}$. This prevents trees from getting overly complex and provides regularization (here, we rely on the default values $\alpha=0.95$, and $\beta=2$, which perform well across many datasets). For the splitting variables, we choose a uniform prior. This implies that each predictor is equally likely to be selected as a splitting variable. We further assign a uniform prior to all thresholds within the splitting rules, based on the range of the respective splitting variable. Next we specify the prior for the terminal node parameters. On these parameters $\bm{\mathrm{m}}_{is,l}$, for $l = 1,\hdots,\#\text{TN}_{is}$, where $\#\text{TN}_{is}$ denotes the number of terminal node parameters of tree $s$ in equation $i$, we impose independent conjugate Gaussian priors that are symmetric across trees and identical for all terminal nodes. As suggested by chipman2010bart the moments of these priors are chosen in a data-driven manner, such that $95$% of the prior probability lies in the interval $(\min(\bm{y}_i),\max(\bm{y}_i))$, where $\bm{y}_i = (y_{i1},\hdots,y_{iT})'$, and such that shrinkage increases the more trees $S$ are chosen for estimation. We choose $S=250$ trees which has been shown to work well for typical macroeconomic time series applications huber2023nowcasting. Additional details are provided in Appendix (ref). \subsection{Other Priors and Sampling Algorithm} Inspired by carriero2021addressing, we assume that $\bm{\Sigma}_t = s_t^2 \bm{\Sigma}$ --- the covariance structure varies proportionally over time and $s_t$ is used to capture outliers. Our prior setup for the constant part of the covariance matrix follows esser2024seemingly. Specifically, we use a hierarchical inverse Wishart prior $\bm{\Sigma} \:\vert\:\{a_i\}_{i=1}^n \sim \mathcal{W}^{-1}(s_0, \bm{S}_0),$ where $s_0 = \nu + n - 1$, $\bm{S}_0 = 2 \nu \cdot \text{diag}(1/a_1,\hdots,1/a_n)$ and $a_i\sim\mathcal{G}^{-1}(1/2,1/A_j^2)$ for $i = 1,\hdots,n,$ and a fixed scale parameter $A_j > 0$. Setting $\nu = 2$ implies a comparatively uninformative prior about the implied correlation structure, different from fixed-hyperparameter versions of this prior which has a tendency to overshrink. If applicable, we assume that $\Pr(s_t = 1) = 1 - \mathfrak{p}$ and $\Pr(s_t \sim \mathcal{U}(2,\overline{\mathfrak{s}})) = \mathfrak{p}$, where $\mathcal{U}(2,\overline{\mathfrak{s}})$ is a discrete uniform distribution with (integer) support between $2$ and $\overline{\mathfrak{s}} = 6$ and $\mathfrak{p}\sim\mathcal{B}(a_\mathfrak{p},b_\mathfrak{p})$ is the probability of observing an outlier. We set $a_\mathfrak{p} = 1$, $b_\mathfrak{p} = 50$ which a priori implies about $2$% of the observations are outliers. Alternative and more flexible stochastic volatility specifications are available in this context chan2023comparing.\footnote{Another computationally attractive variant is to assume $\bm{\Sigma}_t = \bm{S}_t\bm{\Sigma}\bm{S}_t$ where $\bm{S}_t = \text{diag}(s_{1t},\hdots,s_{nt})$, which allows for variable-specific outlier detection using the same setup as in Section (ref). These assumptions allow to invert $\bm{\Sigma}$ only once instead of having to invert $\bm{\Sigma}_t$ for all $t$ to compute the moments in ((ref)).} We obtain draws from the joint posterior of our model using a fairly straightforward MCMC algorithm. Details are provided in Appendix (ref). \section{Empirical Applications} We employ the proposed framework in three related yet distinct applications. First, we use the annual stress test scenarios conducted by the Federal Reserve System and compute forecasts conditional on multiple observables. Second, we study the implications of varying financial conditions on tail risks of output growth, inflation, and employment. Third, we identify a US-based financial shock and gauge the role of spillovers and spillbacks.\footnote{We sometimes compare a \textit{hom}oskedastic (\texttt{BART-hom}, $s_t = 1$ for all $t$) and \textit{het}eroskedastic BART (\texttt{BART-het}) with a \textit{het}eroskedastic BVAR (\texttt{BVAR-het}), featuring the outlier specification. The linear BVAR serves for comparisons as it is a popular workhorse model in related contexts crump2021large; due to our sample featuring the Covid-19 pandemic, we disregard its homoskedastic version.} \subsection{Stress Testing Scenarios for the US Economy} In our first application, we conduct a scenario analysis for the US economy inspired by the 2025 version of the \href{https://www.federalreserve.gov/publications/dodd-frank-act-stress-test-publications.htm}{\textit{Dodd-Frank Act}} (DFA) stress test assumptions. This annual exercise is conducted and published by the Board of Governors of the Federal Reserve System. Details about the underlying dataset are provided in Appendix (ref). The information set features about $25$ broad variables (capturing economic activity, labor market, prices, housing and the financial sector). We estimate our models using quarterly data ranging from 1976Q1 to 2024Q4, and subsequently consider a baseline and adverse scenario for the period from 2025Q1 to 2027Q4. These scenarios are imposed via constraints on the path of the unemployment rate (\texttt{UNRATE}), CPI inflation (\texttt{CPIAUCSL}), and $10$-year government bond yields (\texttt{GS10}), inspired by chan2023conditional. We set the tightness of the restrictions in $\bm{\Omega}_{1:H}$ using the marginal variances on the diagonal of $\bm{\Sigma}$ in each MCMC draw. \begin{figure}[t] \caption{Conditional forecasts with \texttt{BART-het} for selected variables.} \caption*{ \textit{Notes}: Posterior median alongside 50/68 percent credible sets. Restricted variables: consumer price inflation (\texttt{CPIAUCSL}), unemployment rate (\texttt{UNRATE}), $10$-year government bond yields (\texttt{GS10}); unrestricted variables: Real GDP (\texttt{GDPC1}), industrial production (\texttt{INDPRO}), personal consumption expenditure inflation (\texttt{PCECTPI}), payroll employment (\texttt{PAYEMS}), federal funds rate (\texttt{FEDFUNDS}) and excess bond premium (\texttt{EBP}).} \end{figure} The posterior median forecasts and $50$ and $68$ percent posterior credible sets, obtained with \texttt{BART-het} (conditional forecast distributions for other specifications are in Appendix (ref)), for the restricted and selected unrestricted variables, are shown in Figure (ref): real GDP (\texttt{GDPC1}), industrial production (\texttt{INDPRO}), personal consumption expenditure (PCE) inflation (\texttt{PCECTPI}), payroll employment (\texttt{PAYEMS}), federal funds rate (\texttt{FEDFUNDS}) and the gilchrist2012credit excess bond premium (\texttt{EBP}). The baseline scenario draws from the consensus projections from 2025 \textit{Blue Chip Financial Forecasts} and \textit{Blue Chip Economic Indicators}; the adverse scenario is characterized by a recession. Figure (ref) displays corresponding UGIRFs across model specifications, which are computed as the indicated scenario minus the unconditional forecast. The unconditional forecasts from \texttt{BART-het} approximately coincide with the baseline scenario (the UGIRF credible sets cover $0$ in most cases). For the adverse scenario, a different picture emerges --- the scenario forecasts differ significantly from the unconditional forecasts. There is an immediate downturn of economic activity and a reduction in payroll employment. Financial conditions tighten initially but tend to improve subsequently, partially through a monetary easing response by the central bank as reflected in the policy rate. In addition, the assumed trajectory of the conditioning variables results in a modestly disinflationary episode that vanishes by $2027$. \begin{figure}[ht] \begin{subfigure}{\linewidth} \caption{\texttt{BART-het}} \end{subfigure} \begin{subfigure}{\textwidth} \caption{\texttt{BART-hom}} \end{subfigure} \begin{subfigure}{\textwidth} \caption{\texttt{BVAR-het}} \end{subfigure} \caption{Unorthogonalized GIRF (UGIRF, scenario minus unconditional forecast) across model specifications for selected variables.} \caption*{ \textit{Notes}: Posterior median alongside 50/68 percent credible sets; adverse and baseline scenarios in red and blue, respectively. Restricted variables: consumer price inflation (\texttt{CPIAUCSL}), unemployment rate (\texttt{UNRATE}), $10$-year government bond yields (\texttt{GS10}); unrestricted variables: Real GDP (\texttt{GDPC1}), industrial production (\texttt{INDPRO}), personal consumption expenditure inflation (\texttt{PCECTPI}), payroll employment (\texttt{PAYEMS}), federal funds rate (\texttt{FEDFUNDS}) and excess bond premium (\texttt{EBP}).} \end{figure} The conditional forecast distributions across model specifications are similar for most variables apart from output growth, industrial production, and employment. For \texttt{BART-het}, the magnitude of the peak contraction reduces by about half. This can be explained by noting that the algorithm decides to classify several observations (that are otherwise informative about directional movements of variables when assuming homoskedasticity) as outliers. We note that BART, due to the way how tree-based approaches fit data, is capable of dealing with outliers and heteroskedastic data features in the conditional mean function by design huber2023nowcasting,clark2021tail, even when assuming constant variances. But in this case several observations are classified as noise rather than signal. For linear models, the forecast trajectories are smoother, and an overshooting behavior is noticeable around $2027$. Moreover, the response of the federal funds rate is somewhat more persistent, which may be due to nonlinearities arising at the effective lower bound which the BVAR fails to capture. \FloatBarrier \subsection{Financial Conditions in the US and Tail Risk Scenarios} In the next empirical application, we restrict our sample to 1976Q1--2017Q4 and consider the period from 2018Q1 until 2019Q1 as a laboratory to assess nonlinearities between economic variables and financial conditions. For this application we use \texttt{BART-hom} (since the pandemic observations are excluded). We investigate nonlinear patterns of macroeconomic risk, which, following the “growth-at-risk” approach of adrian2019vulnerable, is defined as the predictive quantiles of some variable of interest at a pre-defined probability level (in line with value-at-risk, VaR, in finance). We pick this period because the information set already contains the global financial crisis (and the model thus had the opportunity to learn from this severe financial episode), and because this “holdout sample” otherwise coincides with a comparatively eventless period. We impose hard constraints on the $h = 1$ value of the National Financial Conditions Index (NFCI) and trace the effects of these scenarios on several macroeconomic variables. The scenarios are defined to reflect an increase of the NFCI by approximately $1$, $3$ and $6$ unconditional standard deviations (SDs, reflecting tighter financial conditions) in 2018Q1, i.e., in $\mathcal{C}_1$, which we implement by placing these values as hard restriction on the NFCI in that quarter. From 2018Q2 onward we leave the future unrestricted, i.e., $\mathcal{C}_h = \emptyset$ for $h > 1$. We investigate growth-at-risk (quantiles of real GDP), inflation-at-risk (quantiles of PCE inflation) and labor-at-risk (quantiles of growth in payroll employment) as our objects of interest adams2021forecasting,pfarrhofer2022modeling,clark2024investigating,lopez2024inflation. The CF distributions (density estimates), for average growth rates $\overline{\bm{y}}_{T+h} = \sum_{j=1}^h\bm{y}_{T+j} / h$, are shown in Figure (ref). In the lower panel, we show the difference between the conditional scenario distributions and the unconditional one, which yields the UGIRF (cumulated for all variables except the NFCI). \begin{figure}[t!] \caption{Conditional forecast distributions for selected variables and macroeconomic value-at-risk (VaR) and unorthogonalized GIRFs for different scenarios of financial stress.} \caption*{ \textit{Notes}: “Max. NFCI” refers to the maximum value of the NFCI for the scenarios, with moderate (1), severe (3, comparable to the global financial crisis), and extreme stress (6). Sampling period 1976Q1 to 2017Q4, hard restriction applies 2018Q1. Posterior median alongside 50/68 percent credible set. Variables: Real GDP (\texttt{GDPC1}), personal consumption expenditure inflation (\texttt{PCECTPI}), payroll employment (\texttt{PAYEMS}), national financial conditions index (\texttt{NFCI}). Distribution of average growth rates from 2018Q1 until the quarter indicated on the y-axis in the upper panel. Cumulated UGIRFs for all variables except the NFCI.} \end{figure} The different NFCI scenarios for 2018Q1 shift the predictive distributions. In all cases, the economy contracts which is reflected in a decrease of real GDP growth and payroll employment, and the simulated shock has a modestly disinflationary effect. While the upper tails of the distributions (upside risk) remain comparatively stable, downside risk as measured by the lower quantiles increases significantly for all considered variables (the red shaded VaR $<0.05$ moves strongly leftwards), and there are some visible asymmetries. While the moderate NFCI scenario (max. NFCI $=1$) results in growth-at-risk for the $5$th percentile at about $-5$ percent, the severe (max. NFCI $=3$) and extreme (max. NFCI $=6$) stress scenarios yield $-7$ and $-12.5$ percent, respectively. This finding is also present for inflation-at-risk and labor-at-risk. The resulting predictive distributions exhibit non-Gaussian features, chief among them being heavy tails and skewness clark2021tail. In addition, there are hints of multimodality as the assumed values for the NFCI in 2018Q1 turn more extreme, which relates to the discussions in adrian2021multimodality. These features can arise --- even in one-step ahead predictions and for a single restricted period --- due to, for example, the initial conditions of the nonlinear unconditional mean function at the forecast origin. \FloatBarrier \subsection{Spillovers and Spillbacks of US Financial Shocks} In our final application, we estimate the effect of a financial shock in the US and trace its effects through the domestic economy, but also capture spillovers and spillbacks to and from other economies. We use an adjusted dataset in this case, which drops several of the domestic indicators, but adds bilateral exchange rates alongside real GDP for the EA and the UK. To identify the financial shock, we place timing-restrictions on the contemporaneous impulse responses. This is operationalized with a specific ordering of the quantities in the vector $\bm{y}_t$ --- we structure this vector such that all slow moving domestic and foreign macroeconomic variables come first (which imposes zero restrictions on impact). These variables are then followed by the EBP, and all fast moving variables such as those capturing the financial economy. We then use a Cholesky decomposition of the form $\bm{\Sigma} = \bm{P}\bm{P}'$ where $\bm{P}$ is lower triangular. That is, $\bm{H}^{-1} = \bm{P}$, and we orthogonalize the structural shocks different to our previous applications. The orthogonalized innovation of the EBP equation is interpreted as the financial shock, similar to gilchrist2012credit,barnichon2022effects. In a multicountry context, huber2024asymmetries use an identification scheme identical to ours. We use $d_0 = 0$ and simulate different shock sizes and signs with $d \in \{-3,-1,1,3,6\}$. In contrast with conventional linear frameworks, our approach allows to assess nonlinearities of higher-order responses with respect to different signs and magnitudes of a proportional shock impact with GIRFs. Such asymmetries and related nonlinearities have recently gained attention both in a VAR and local projection context, see, e.g., mumtaz2022impulse,carriero2023shadow,forni2024nonlinear,hauzenberger2024machine. While our framework allows to compute dynamic responses for each period in our sample, we focus on time averages in the results that follow. Further, we rescale all SGIRFs by computing $\bm{\delta}_{\tau,h}^{(d)} / d$ so that all responses shown below reflect a $1$ SD financial shock gonccalves2024state,kolesar2024dynamic. Note that for linear VARs, such scaling yields identical IRFs for all shock sizes. For nonlinear models, this is not necessarily the case and allows for a visual inspection of asymmetries. Selected variables are shown in Figure (ref). The rows in the figure show different subsets of the same results, structured such that the shocks of different signs and sizes can be compared with ease. \begin{figure}[!ht] \caption{Structural GIRFs for selected variables to a financial shock in the US, comparing asymmetries due to size and sign of the shocks.} \caption*{ \textit{Notes}: Posterior median alongside 50/68 percent credible sets. Cumulated responses for variables in differences and levels for all other variables. “Horizon” refers to periods after impact of the shock. Variables: Real GDP (\texttt{GDPC1}), consumer price inflation (\texttt{CPIAUCSL}), payroll employment (\texttt{PAYEMS}), federal funds rate (\texttt{FEDFUNDS}) and S&P500 index (\texttt{SP500}).} \end{figure} Starting with the first row, we find that the size of the financial shock causes limited asymmetries in responses for the indicated variables. Financial shocks of different sizes rather symmetrically decrease real GDP and payroll employment. Interestingly, the effects of very large shocks increase slightly less than proportionally. Peak effects occur about two years after impact of the shock. In addition, the shock puts a persistent downward pressure on prices and leads to a decline in the federal funds rate which peaks at about $-15$ basis points after around a year. Notably the federal funds rate does not react on impact, different to stock returns which immediately decline by about $1.5$ percent during the quarter when the shock materializes. Qualitatively and in terms of magnitudes, these estimates are roughly in line with the previous literature. The US-based financial shock spills over to the other economies and leads to contractionary effects in terms of real GDP, see Appendix (ref). Having established that the size of the financial shock does not seem to matter much, the second and third row zoom into sign asymmetries. The $1$ SD US-based financial shock does not result in significant sign asymmetries (the posterior distributions overlap for the most part). Turning to the final row, this clearly differs for larger sized shocks of different signs. For these GIRFs that show the responses to a positive (adverse) and negative (benign) $3$ SD shock, asymmetries are visible for most variables. In particular, we find that the negative effect on payroll employment is almost twice as large for adverse shocks, and the Federal Reserve responds more strongly to adverse financial shocks, as measured by the much stronger shift in the federal funds rate. These findings corroborate previous evidence forni2024nonlinear,hauzenberger2024machine. \begin{figure}[t] \caption{Structural and restricted GIRFs for selected variables.} \caption*{ \textit{Notes}: The restricted case assumes that the US financial shock does not spill over to non-domestic variables. Posterior medians alongside 50/68 percent credible sets. Cumulated responses for variables in differences and levels for all other variables. “Horizon” refers to periods after impact of the shock. Variables: Real GDP (\texttt{GDPC1}), consumer price inflation (\texttt{CPIAUCSL}), payroll employment (\texttt{PAYEMS}), federal funds rate (\texttt{FEDFUNDS}) and S&P500 index (\texttt{SP500}).} \end{figure} We next explore the role of international variables in the domestic transmission of the US shock. For this purpose, besides SGIRFs, we consider an alternative analysis where non-domestic transmission channels are switched off --- we investigate how international channels affect the domestic transmission of shocks originating in the US. That is, we impose the restriction that foreign variables (real GDP in the EA and UK, and exchange rates) do not respond to the financial shock in the US in this counterfactual, thereby simulating a scenario where the financial shock is confined to the domestic economy without any real or financial spillovers (or spillbacks); see breitenlechner2022goes for a monetary application in this context. We restrict all domestic shocks to their unconditional distribution, and use the non-domestic ones as driving shocks to compute the RGIRFs. The results are displayed in Figure (ref). The upper panels show the SGIRFs (those shown and discussed in the context of Figure (ref)) and “no spillovers” RGIRFs. Two key findings are worth reporting. First, for the most part, ruling out spillovers and spillbacks does not significantly alter the dynamic responses after an adverse shock. It is worth mentioning, however, that restricting the international transmission leads to slightly smaller effects on average. Second, international transmission channels appear to matter for inflation dynamics, and especially so for benign financial shocks. The response of inflation turns insignificant in this case, which is also associated with a less forceful action by the central bank as captured in more muted response of the policy rate. \FloatBarrier \section{Conclusions} This paper presents a unified methodology for conducting scenario analysis in multivariate macroeconomic settings, accommodating nonlinearities and unknown functional forms of conditional mean relationships. These methods are applicable to traditional nonlinear frameworks, such as variants of threshold or time-varying parameter models, but also to more recently developed models incorporating Bayesian machine learning. Our framework addresses some limitations of linear and parametric models in generating various types of counterfactual analyses and is suitable for large macroeconomic datasets. The empirical applications, using Bayesian additive regression trees as an example of nonparametric modeling of the conditional mean function, underscore the role of nonlinearities in shaping macroeconomic dynamics. For instance, a scenario analysis based on Federal Reserve stress test assumptions reveals differences between linear and nonlinear models in forecasting economic contractions and recoveries. Similarly, in a growth-at-risk application we measure nonlinear macroeconomic risks under financial stress. Finally, an analysis of financial spillovers reveals asymmetries in the transmission of shocks and the influence of international linkages on domestic outcomes. {\setstretch{1.2}\putbib}