Bayesian Forecasting in Economics and Finance: A Modern Review

As a general rule, if $p(\bm{\theta}\mathbf{|y}_{1:T})$ does not have a closed-form representation, it is also not amenable to Monte Carlo sampling, as the latter requires that $p(\bm{\theta}\mathbf{|y}_{1:T})$ can be decomposed into recognizable densities, from which computer simulation is feasible. IS (Kloek and van Dijk,, 1978; Geweke,, 1989), via use of an ‘importance’ or ‘proposal’ density, $q(\bm{\theta}\mathbf{|y}_{1:T})$, that matches $p(\bm{\theta}\mathbf{|y}_{1:T})$ well and which can be drawn from, is a possible solution in some cases. However, the algorithm can fail to produce representative draws from $p(\bm{\theta}\mathbf{|y}_{1:T})$ when the dimension of $\bm{\theta}$ is large, due to the difficulty of finding a $q(\bm{\theta}\mathbf{|y}_{1:T})$ that is a ‘good match’ to $p(\bm{\theta}\mathbf{|y}_{1:T})$ in high dimensions.

In contrast, under certain conditions, a Gibbs sampler is able to produce a (dependent) set of draws from the joint posterior via iterative sampling from lower dimensional, and often standard, conditional posteriors. In other words, a Gibbs sampler takes advantage of the fact that, while joint and marginal posterior distributions are usually complex in form and unable to be simulated from directly, conditional posteriors are often standard and amenable to simulation. Given the satisfaction of the required convergence conditions (Geyer,, 2011), draws $\bm{\theta}^{(i)}$, $i=1,2,...,M$, produced via iterative sampling from the full conditionals converge in distribution to $p(\bm{\theta}\mathbf{|y}_{1:T})$ as $M\rightarrow\infty$, and can be used to produce a $\sqrt{M}$-consistent estimate of the ordinates of $p(y_{T+1}|\mathbf{y}_{1:T})$ across the support of $y_{T+1}$ in the manner described in Section 2.2. Decisions about how to partition, or ‘block’ $\bm{\theta}$ need to be made (Liu et al.,, 1994; Roberts and Sahu,, 1997), with a view to increasing the ‘efficiency’ of the chain which, in effect, amounts to ensuring an accurate estimate of $p(y_{T+1}|\mathbf{y}_{1:T})$ for a given number of draws, $M.$ (See Algorithm 1 in Appendix A.1.)

Chib, (1993) and McCulloch and Tsay, (1994) are the earliest examples of using Gibbs algorithms for Bayesian estimation and prediction in time series settings; both papers exploiting the fact that despite $p(\bm{\theta}\mathbf{|y}_{1:T})$ and $p(y_{T+1}|\mathbf{y}_{1:T})$ precluding analytical treatment in most of the examples considered, the conditional posteriors always have closed forms. As one would anticipate however, a ‘pure’ Gibbs algorithm based on a full set of standard conditionals is not always possible, with the more typical situation being one in which one or more of the conditionals – associated with any given partitioning of the parameter space – are not available in closed form. The following section describes how to adapt a Gibbs algorithm in cases where certain conditional components are not known in closed form, and, in so doing, illustrates a powerful simulation-based algorithm for accessing $p(y_{T+1}|\mathbf{y}_{1:T})$ in more complex settings.

The Gibbs sampler is only one example of an MCMC algorithm. The first such example – the ‘Metropolis’ algorithm – appeared in a paper that has assumed an important status in the history of statistics: Metropolis et al., (1953)³³3For example, Dongarra and Sullivan, (2000) rank the Metropolis algorithm proposed in Metropolis et al., as one of the 10 algorithms “with the greatest influence on the development and practice of science and engineering in the 20th century”.. The Metropolis algorithm was subsequently generalized by Hastings, (1970), and it is this ‘MH’ version of the method that is typically referenced. For the purpose of this review, the key role of the MH algorithm is to enable sampling from non-standard conditionals within a Gibbs algorithm, in particular when the dimension of the conditionals precludes (say) the exclusive use of inverse cumulative distribution function (ICDF) sampling.⁴⁴4Any non-standard probability distribution can, in principle, be drawn from using ICDF sampling. The term ‘Griddy Gibbs’ sampling was first used by Ritter and Tanner, (1992) to refer to the use of ICDF sampling to draw from non-standard conditionals in a Gibbs scheme. Given that the method amounts to the use of numerical quadrature, it suffers from the curse of dimensionality, and is thus infeasible for drawing from anything other than very low-dimensional conditionals. See Bauwens and Lubrano, (1998) for the application of the Griddy-Gibbs sampler to a generalized autoregressive conditionally heteroscedastic (GARCH) model for financial returns.

Under regularity, a Markov chain that converges to $p(\bm{\theta}\mathbf{|y}_{1:T})$ can be produced by embedding an MH algorithm (or MH algorithms) within an outer Gibbs loop. In short, an MH-within-Gibbs algorithm proceeds by drawing from any non-standard conditional indirectly, via a ‘candidate’, or ‘proposal’ distribution that is deemed to be a good match to the inaccessible conditional, and accepting the draw with a given probability. Critically, the formula that defines the acceptance probability involves evaluation of the non-standard conditional only up to its integrating constant; hence the conditional need not be known in its entirety.⁵⁵5Moreover, and in contrast to IS, the requirement to find a well-matched proposal distribution is facilitated by the dimension reduction invoked by the breaking down of the high-dimensional joint posterior into the lower dimensional conditionals, before any proposal distribution needs to be specified. Again, under appropriate regularity, the draws $\bm{\theta}^{(i)}$, $i=1,2,...,M$, from the MH-within-Gibbs algorithm converge in distribution to $p(\bm{\theta}\mathbf{|y}_{1:T})$ as $M\rightarrow\infty$, and can be used to produce a $\sqrt{M}$-consistent estimate of the ordinates of $p(y_{T+1}|\mathbf{y}_{1:T})$. (See Algorithm 2 in Appendix A.2.)

As will become evident in the subsequent empirical review sections, MH-within-Gibbs algorithms remain the dominant form of method used to sample from posteriors – and to estimate predictive distributions – for time series models for which a convenient partitioning of the parameter space is available, and for which the conditional posteriors are known up to their integrating constants. Hence, we reserve further elaboration on the use of such algorithms in practice until the appropriate points in Section 4.

For many empirical problems in economics and related fields, a suitable model can be partitioned into two sets: static unknowns $\bm{\theta}$, which are fixed throughout time, and latent data, $\mathbf{z}_{1:T}=(z_{1},z_{2},...,z_{T})^{\prime}$, which vary over time. The latent states may be intrinsic to the model – as in a state space model – or may be auxiliary variables introduced purely for the purpose of facilitating posterior sampling. Application of a Gibbs-based MCMC scheme to the joint, or ‘augmented’ set of unknowns $\left(\bm{\theta},\mathbf{z}_{1:T}\right)$ is often referred to as ‘data augmentation’, in the spirit of Tanner and Wong, (1987), and such schemes have enabled the Bayesian analysis of large classes of time series models that would otherwise have been inaccessible.

We illustrate here the basic principles of the approach using a state space model governed by a measurement density for the observed scalar random variable, $y_{t}$, and a Markov transition density for a scalar state variable, $z_{t}$,

Using the generic notation in (7) and (8), the augmented posterior is

In certain cases, the model structure is such that a pure Gibbs scheme can be used to produce draws from $p(\bm{\theta},\mathbf{z}_{1:T}|\mathbf{y}_{1:T})$ and, thus, from $p(\bm{\theta}|\mathbf{y}_{1:T})$; an insight obtained independently by Carter and Kohn, (1994) and Frühwirth-Schnatter, (1994) for the case of the linear Gaussian state-space model, for example. However, implementation of such a scheme will, by definition, require both $p(\bm{\theta}|\mathbf{z}_{1:T},\mathbf{y}_{1:T})$ and $p(\mathbf{z}_{1:T}|\bm{\theta},\mathbf{y}_{1:T})$ to have recognizable forms. In more general cases, in which either the measurement or state equation has non-linear and/or non-Gaussian features, the resulting conditionals will not necessarily have a known closed form, which necessitates the addition of MH steps within the outer Gibbs loop. Such a treatment was the method of attack for large classes of models in the 1990s and early 2000s. Relevant contributions here, which include specific treatments of the ubiquitous stochastic volatility (SV) model, are Polson et al., (1992), Jacquier et al., (1994), Shephard and Pitt, (1997), Kim et al., (1998), Chib et al., (2002), Stroud et al., (2003), Chib et al., (2006), Strickland et al., (2006), Omori et al., (2007) and Strickland et al., (2008). The reviews of Fearnhead, (2011) and Giordani et al., (2011) provide more detailed accounts and extensive referencing of this earlier literature.⁶⁶6We also note here the work of Chib and Greenberg, (1994), in which the state space representation of an autoregressive moving average (ARMA(p,q)) model (Harvey,, 1981) was exploited, and the principle of data augmentation invoked, in order to enable an MH-within-Gibbs scheme to be applied. (See also Appendix A.3.)

To conclude, and once again using the generic notation in (7) and (8), once draws have been produced from $p(\bm{\theta},\mathbf{z}_{1:T}|\mathbf{y}_{1:T})$, the predictive pdf,

can be estimated in the usual way, using subsequent draws from $p(z_{T+1}|z_{T},\bm{\theta})$ and $p(y_{T+1}|z_{T+1},\bm{\theta}\mathbf{,y}_{1:T})$, or by averaging the conditional predictives over all draws of $z_{T+1}$ and $\bm{\theta}$.

The MCMC methods that evolved during the late 20th century continue to serve as the ‘bread and butter’ of Bayesian forecasting, as will be made evident in Section 4. Nevertheless, more ambitious forecasting problems are now being tackled, and this has tested the mettle of some of the early algorithms. As a consequence, Bayesian forecasters have begun to exploit more modern computational techniques, and it is those techniques that we touch on briefly in this section.

It is convenient to characterize these newer computational developments as different types of solutions to so-called ‘intractable’ forecasting problems, by which we mean: (a) forecasts based on models with data generating processes (DGPs) that cannot be readily expressed as a pdf, or probability mass function (pmf); (b) forecasts based on high-dimensional models, with a very large number of unknowns; (c) forecasts produced using extremely large data sets. Problems that feature problem (a) are referred to as doubly-intractable problems, as not only is $p(\bm{\theta}|\mathbf{y}_{1:T})$ not available in its entirety (as is typical), but the DGP itself is also not able to be expressed analytically.

With reference to (a), the MCMC methods referenced so far entail the evaluation of the DGP as a pd(/m)f, either in the calculation of the acceptance probability in any MH sub-step, or in the specification of full conditionals in any ‘pure’ Gibbs step. Hence, they are infeasible when DGPs do not admit such a representation. Many such DGPs exist (see, for example, Martin et al., 2023a, , for a list of examples); however, particularly pertinent ones to mention here are continuous time models in finance with unknown transition densities (Gallant and Tauchen,, 1996), $\alpha$-stable models for financial returns (and/or their volatility) (Peters et al.,, 2012; Martin et al.,, 2019), and stochastic dynamic equilibrium models in economics (Calvet and Czellar,, 2015). With regard to (b), whilst, in principle (and under appropriate regularity), a convergent MCMC chain can be constructed for any model, the exploration of a very high-dimensional parameter space via an MCMC algorithm can be prohibitively slow (Tavaré et al.,, 1997; Rue et al.,, 2009; Braun and McAuliffe,, 2010; Lintusaari et al.,, 2017; Betancourt,, 2018; Johndrow et al.,, 2019). Hence, in models with a very large number of unknowns – including those with multiple sets of latent variables – the production of an accurate MCMC-based estimate of $p(y_{T+1}|\mathbf{y}_{1:T})$ in a practical amount of time may not be possible. Finally, regarding (c) MCMC schemes require pointwise (i.e. for each $y_{t}$) evaluation of $p(\mathbf{y}_{1:T}|\bm{\theta})$ at each draw of $\bm{\theta}$, thereby inducing an $O(T)$ computational burden at each iteration in an MCMC chain.⁷⁷7We recall that a sequence $X_{T}$ is $O(T)$ if $|X_{T}/T|$ is bounded as $T\rightarrow+\infty$. Such schemes can thus struggle when confronted with ‘big data’ (Bardenet et al.,, 2017). In this context, ‘big data’ refers to situations where, due to the length and/or size of the data set, the repeated evaluation of the likelihood function that is required to produce draws from the corresponding MCMC chain is too time consuming for the algorithm to run in a reasonable amount of time.

The methods in the following sections have been designed to solve one or more of these instances of intractability. The techniques in Section 3.2.2 do so whilst preserving the ‘exact’ nature of the estimate of $p(y_{T+1}|\mathbf{y}_{1:T})$, whilst those in Section 3.2.3 aim to produce an approximation of $p(y_{T+1}|\mathbf{y}_{1:T})$ only.

The first two decades of the 21st century have witnessed a wealth of advances in both MCMC and IS-based algorithms. The goal of the newer MCMC algorithms – at their heart – is to explore the high mass region of the joint posterior more efficiently, in particular when the dimension of the space of unknowns is large. This, in turn, enables a more accurate estimate of $p(y_{T+1}|\mathbf{y}_{1:T})$ to be produced for a given computational budget. This goal has been achieved via a variety of means, which (in the spirit of Robert et al.,, 2018, and Martin et al., 2023b, ) can be summarized as follows: the use of more geometric information about the target posterior, most notably the use of Hamiltonian updates (Neal, 2011b, ; Hoffman and Gelman,, 2014); the use of better MH candidate, or proposal distributions, including those that ‘adapt’ to previous draws (Nott and Kohn,, 2005; Roberts and Rosenthal,, 2009); various types of combinations of multiple chains (Jacob et al.,, 2011; Neal, 2011a, ; Neiswanger et al.,, 2013; Glynn and Rhee,, 2014; Huber,, 2016; Jacob et al.,, 2020); or the use of ex-post variance reduction methods (Craiu and Meng,, 2005; Douc and Robert,, 2011; Owen,, 2017; Baker et al.,, 2019). We refer the reader to Green et al., (2015), Robert et al., (2018) and Dunson and Johndrow, (2020) for detailed reviews of modern developments in MCMC, and to Jahan et al., (2020) for an overview of the way in which certain of the newer methods manage the problem of scale – in terms of either the unknowns or the data, or both.

Whilst not designed expressly to deal with problems of scale, sequential Monte Carlo (SMC) methods – which exploit the principles of IS – have developed in parallel to the expansion of the MCMC stable. Devised initially for the sequential analysis of state space models, via methods of ‘particle filtering’ (Gordon et al.,, 1993), SMC methods have evolved into a larger suite of methods used to perform both sequential and non-sequential tasks (Naesseth et al.,, 2019; Chopin and Papaspiliopoulos,, 2020). For the purpose of this review, the most pertinent development is the melding of particle filtering with MCMC in state space settings to produce a particle marginal MH (PMMH) algorithm (Andrieu et al.,, 2011; Flury and Shephard,, 2011; Pitt et al.,, 2012; Doucet et al.,, 2015; Deligiannidis et al.,, 2018). Such algorithms tackle intractability type (a) in the dichotomy of the previous section, by replacing an ‘unavailable’ likelihood function by an unbiased estimate – produced via the particle filter – in an MH algorithm which, under regularity, retains the posterior $p(\bm{\theta}|\mathbf{y}_{1:T})$ as its invariant distribution. Given the increasingly important role played by PMMH, a brief algorithmic description of it is included in Algorithm 3 in Appendix A.4.⁸⁸8PMMH is actually a special case of the general pseudo-marginal MH technique (also sometimes denoted by the abbreviation ‘PMMH’), in which a pseudo likelihood, produced – in some manner or another – as an unbiased estimator of the true likelihood, is used within an MH algorithm. See, for example, the subsampling methods based on pseudo-marginal MCMC (Bardenet et al.,, 2017; Quiroz et al.,, 2018; Quiroz et al.,, 2019) used expressly to improve the performance of MCMC in the case of a large-dimensional $\mathbf{y}_{1:T}$ (i.e. intractability type (c)).

In situations in which the dimension, or structure of the forecasting model, or the size of the data set, still precludes the use of either an MCMC or a PMMH approach, an approximate method may be the only computational option. The cost of adopting such a solution is that these methods no longer directly target the exact predictive, $p(y_{T+1}|\mathbf{y}_{1:T})$; instead, an approximation of $p(y_{T+1}|\mathbf{y}_{1:T})$ becomes the goal.

The spirit of these methods is to approximate $p(y_{T+1}|\mathbf{y}_{1:T})$ via some feasible approximation to the posterior $p(\bm{\theta}|\mathbf{y}_{1:T})$. Denoting the posterior approximation generically by $g(\bm{\theta}|\mathbf{y}_{1:T})$, the resultant approximate predictive can be expressed as

in the case where there are only static unknowns. When the model features both static parameters and time-varying latent parameters, and exploiting the Markov property of the state process in (8), the approximate predictive can be represented as

Given draws of $\bm{\theta}$ from $g(\bm{\theta}|\mathbf{y}_{1:T})$, and given an appropriate forward-filtering algorithm to draw from $p(z_{T}|\bm{\theta},\mathbf{y}_{1:T})$ when needed, a simulation-based estimate of $g(y_{T+1}|\mathbf{y}_{1:T})$ can be produced in the usual way, either as a sample mean of the conditional predictives defined by the draws of $\bm{\theta}$ (and $z_{T+1}$), or by applying kernel density techniques to the draws of $y_{T+1}$ from the conditional predictive.

With reference to the taxonomy of intractable problems delineated in Section 3.2, the different methods of producing $g(\bm{\theta}|\mathbf{y}_{1:T})$ (and, hence, $g(y_{T+1}|\mathbf{y}_{1:T})$) can be categorized according to whether they are being used to obviate (a) or to tackle a problem of scale: (b) and/or (c). Both ABC and BSL avoid the need to evaluate the DGP and, hence, are feasible methods in the doubly-intractable settings of category (a). In brief, both methods require only simulation, not evaluation, of the DGP. The approximation of $p(\bm{\theta}|\mathbf{y}_{1:T})$ arises, primarily, from the fact that both methods – in different ways – degrade the information in the full data set, $\mathbf{y}_{1:T}$, to the information contained in a set of summary statistics, $\eta(\mathbf{y}_{1:T})$. As such, the target becomes the so-called ‘partial’ posterior for $\bm{\theta}$, which conditions on $\eta(\mathbf{y}_{1:T})$, rather than $\mathbf{y}_{1:T}$. The quality of the approximation is thus dependent on the informativeness of the summaries, as well as on other forms of approximation invoked in the implementation of the methods. Vanilla versions of both algorithms are provided in Algorithms 4 (Appendix A.5) and 5 (Appendix A.6) respectively.

In contrast to ABC and BSL, VB and INLA still target the exact posterior $p(\bm{\theta}|\mathbf{y}_{1:T})$, but provide approximations that can be computationally convenient when the scale of the empirical problem is large in some sense (so problem (b) and/or problem (c)), often as a consequence of the specification of a high number of latent, or ‘local’, parameters in the model, in addition to the (usually) smaller set of ‘global’ parameters ($\bm{\theta}$ in our notation). Adopting the technique of the calculus of variations, VB produces an approximation of $p(\bm{\theta}\mathbf{|y}_{1:T})$ that is ‘closest’ to $p(\bm{\theta}\mathbf{|y}_{1:T})$ within a chosen variational family, whilst INLA applies a series of nested Laplace approximations (Laplace,, 1774; Tierney and Kadane,, 1986; Tierney et al.,, 1989) to a high-dimensional latent Gaussian model to produce an approximation of $p(\bm{\theta}|\mathbf{y}_{1:T})$. Both VB and INLA exploit state-of-the-art optimization techniques, for the purpose of minimizing the ‘distance’ between $p(\bm{\theta}\mathbf{|y}_{1:T})$ and the variational approximation in the case of VB, and for the purpose of producing the mode of the high-dimensional vector of latent states in the case of INLA. The basic principles of VB and INLA are provided in Appendices A.7 and A.8 respectively.

We refer the interested reader to Martin et al., 2023a for an extensive review of all of these approximate Bayesian methods, as well as more complete coverage of the existing literature, including references to in-depth reviews of specific methods. Martin et al., 2023a also includes discussion of ‘hybrid’ methods that mix and match features of more than one computational technique, with the aim of tackling multiple instances of ‘intractability’ simultaneously.

Regardless of which approximation method is used, the hope is that the resulting approximate predictive $g(y_{T+1}|\mathbf{y}_{1:T})$ performs well relative to the inaccessible exact predictive, and that issue is addressed in certain work cited in the empirical reviews in Section 4.

Inherent in the conventional Bayesian approach to forecasting is the assumption that the process that has generated the observed data tallies with the particular model that underpins the likelihood function. Bayesian model averaging (BMA) – and the resultant predictive in (3) – has evolved as a principled way of catering for uncertainty about the predictive model, and BMA remains a very important technique in the Bayesian toolbox. Nevertheless, underpinning BMA is still the assumption that the true process is spanned by the set of models over which one averages – i.e. that the so-called $\mathcal{M}$-closed view of the world (Bernardo and Smith,, 1994) prevails.

In response to these perceived limitations of the conventional approach, attention has recently been given to producing predictions that are ‘fit for purpose’, by focusing the Bayesian machinery on the specific goals of the predictive analysis at hand. In the following sections we briefly summarize three such approaches, all of which move beyond the conventional likelihood-based Bayesian update, and $\mathcal{M}$-closed paradigm: seeking to produce accurate predictions without recourse to the assumption of correct model specification.

Loaiza-Maya et al., (2021) propose an approach to Bayesian prediction expressly designed for the context of misspecification. In brief, rather than a correct predictive model being assumed, a prior is placed over a class of plausible predictive models. The prior is then updated to a posterior via a sample criterion function that is constructed using a scoring rule (Gneiting and Raftery,, 2007) that rewards the type of predictive accuracy (e.g. accurate prediction of extreme values) that is important for the particular empirical problem being tackled. With a criterion function that explicitly captures predictive accuracy replacing the likelihood function in the Bayesian update, the explicit need for correct model specification is avoided.

Following Gneiting and Raftery, (2007), and using generic notation, for $\mathcal{P}$ a convex class of predictive distributions on $(\Omega,\mathcal{F})$, the predictive accuracy of $P\in\mathcal{P}$ can be assessed using a scoring rule $S:\mathcal{P}\times\Omega\rightarrow\mathbb{R}$. If the value $y$ eventuates, then the positively-oriented ‘score’ of the predictive $P$, is $S(P,y).$ The expected score under the true unknown predictive $P_{0}$ is defined as

A scoring rule is said to be proper relative to $\mathcal{P}$ if, for all $P,G\in\mathcal{P}$, $S(G,G)\geq S(P,G),$ and is strictly proper, relative to $\mathcal{P}$, if $S(G,G)=S(P,G)\iff P=G$. Scoring rules are important mechanisms as they elicit truth telling within the forecasting exercise: if the true predictive $P_{0}$ were known, then in terms of forecasting accuracy as measured by the scoring rule $S(\cdot,\cdot)$ it would be optimal to use $P_{0}$.

Different scoring rules rewards different forms of predictive accuracy (see Gneiting and Raftery,, 2007, Opschoor et al.,, 2017, and Martin et al.,, 2022 for expositions); hence the motivation to drive the update by the score that ‘matters’. Since $P_{0}$ and the expected score $\mathbb{S}(\cdot,P_{0})$ are unattainable in practice, an estimate based on $\mathbf{y}_{1:T}$ is used to define the sample criterion, $S_{T}(\bm{\theta}):=\sum_{t=0}^{T-1}S[p(y_{t+1}|\bm{\theta},\mathbf{y}_{1:t}),y_{t+1}]$, where $p(y_{t+1}|\bm{\theta},\mathbf{y}_{1:t})$ is the pdf associated with a given $P.$ Adopting the exponential updating rule proposed by Bissiri et al., (2016) (see also Giummolè et al.,, 2017, Holmes and Walker,, 2017, Guedj,, 2019, Lyddon et al.,, 2019, and Syring and Martin,, 2019), Loaiza-Maya et al., (2021) define the generalized (or Gibbs) posterior:

for some learning rate $\omega\geq 0$, calibrated in a preliminary step. This posterior explicitly places high weight on – or focuses on – values of $\bm{\theta}$ that yield high predictive accuracy in the scoring rule $S(\cdot,\cdot)$. As such, the process of building a Bayesian predictive as:

is termed ‘focused Bayesian prediction’ (FBP) by the authors. By construction, when the predictive model, $p(y_{T+1}|\bm{\theta}\mathbf{,y}_{1:T})$, is misspecified, (15) will – out-of-sample – often outperform, in the chosen rule $S(\cdot,\cdot)$, the likelihood (or log-score)-based predictive in (2), and this is demonstrated in Loaiza-Maya et al., both theoretically and in extensive numerical illustrations.

Since a positively-oriented score can, equivalently, be viewed as the negative of a measure of predictive loss, FBP can also be referred to as ‘loss-based’ prediction. Such terminology is indeed adopted in Frazier et al., 2022b , in which the principles delineated here are extended to high-dimensional models, and approximations to both $\pi_{w}(\bm{\theta}|\mathbf{y}_{1:T})$ and $p_{w}(y_{T+1}|\mathbf{y}_{1:T})$ based on VB proposed, and validated. We note that the term ‘loss’ as it is used in Loaiza-Maya et al., (2021) and Frazier et al., 2022b refers specifically to predictive loss as quantified by a proper scoring rule. For the application of loss-based Bayesian inference, in which more general forms of loss functions may drive the Bayesian update, we refer the reader to certain of the other literature cited above, namely Bissiri et al., (2016), Holmes and Walker, (2017), Lyddon et al., (2019) and Syring and Martin, (2019).

The predictive distributions within the ‘plausible class’ referenced above may characterize a single dynamic structure depending on a vector of unknown parameters, $\bm{\theta}$, or may constitute weighted combinations of predictives from distinct models, in which case $\bm{\theta}$ comprises both the model-specific parameters and the weights. As such, FBP provides a coherent Bayesian method for estimating weighted combinations of predictives via predictive accuracy criteria, and without the need to assume that the true model is spanned by the set of constituent predictives – an assumption that underpins BMA, as we have noted.

A similar motivation underlies other contributions to the extensive Bayesian literature on estimating combinations of predictives that has now developed (and which rivals the large frequentist literature on forecast combinations that has also evolved⁹⁹9See Hall and Mitchell, (2007), Ranjan and Gneiting, (2010), Geweke and Amisano, (2011) and Gneiting and Ranjan, (2013) for early contributions to the frequentist forecast combination literature, and Wang et al., (2022) for a recent review. We note that whilst Geweke and Amisano, (2011) is not explicitly Bayesian, in terms of estimating the optimal predictive combination, it provides important insights into the connection between the ‘optimal linear pool’ and BMA, and also uses Bayesian numerical methods in the production of some of the constituent forecast distributions.), with predictive performance – quantified by a range of user-specified measures of predictive accuracy – driving the posterior updating of the weights. Indeed, the Bayesian literature, having access as it does to powerful computational tools, has been able to invoke more complex weighting schemes than can be tackled via frequentist (optimization) methods. Notable contributions, including some also driven by the criterion of predictive calibration (Dawid,, 1982; Dawid,, 1985; Gneiting et al.,, 2007), include Billio et al., (2013), Casarin et al., 2015a , Casarin et al., 2015b , Casarin et al., (2016), Pettenuzzo and Ravazzolo, (2016), Aastveit et al., (2018), Bassetti et al., (2018), Baştürk et al., (2019) and Casarin et al., (2023). Once again adopting the language of Bernardo and Smith, (1994), this literature seeks to move Bayesian predictive combinations beyond the $\mathcal{M}$-closed world of BMA to the $\mathcal{M}$-open world that accords with the reality of misspecification.

We complete this section by also highlighting one particular generalization of BMA that aims, not so much to cater for the $\mathcal{M}$-open world but, rather, to remove the fixed-weight restriction that is inherent to BMA. Certain of the references cited above either explicitly allow for the weights attached to the constituent forecasts to evolve over the time period, $t=1,2,....,T$, on which the predictive distribution for $y_{T+1}$ conditions (e.g. Billio et al.,, 2013, and Casarin et al.,, 2023) or implicitly allow for such a possibility (e.g. Loaiza-Maya et al.,, 2021). However, so-called dynamic model averaging (DMA) accommodates time-varying weights via a more direct generalization of BMA, and nests BMA when appropriate settings are activated (see Koop and Korobilis,, 2012, page 875 for an illustration of this). We refer the reader to Raftery et al., (2010) for the initial proposal of DMA, Koop and Korobilis, (2012) for the application of the method to forecasting inflation, and Nonejad, (2021) for a recent review of the methodology, with a focus on applications in economics and finance.¹⁰¹⁰10We also refer the reader to Green, (1995), Madigan and Raftery,, 1995, and George, (2000), for alternative approaches to catering for model uncertainty in the Bayesian framework, In brief, such approaches – in one way or another – design MCMC samplers to tackle an augmented space in which model uncertainty is incorporated. As a consequence, the computation of any expectation of interest, including that which defines a predictive distribution, automatically factors in all uncertainty associated with both the parameters of each model and the model structure itself. See Green, (2003), Marin et al., (2005), Chib, (2011) and Fan and Sisson, (2011) for reviews and more complete referencing.

A third approach that seeks to produce Bayesian predictions without relying explicitly on correct model specification is Bayesian predictive synthesis (BPS) (Johnson,, 2017; McAlinn and West,, 2019; McAlinn et al.,, 2020; Aastveit et al.,, 2023), recently expanded to Bayesian predictive decision synthesis (BPDS) by Tallman and West, (2022). In particular, BPDS provides a sound decision-theoretic framework for constructing forecast combinations, and can be shown to encompass several commonly-suggested Bayesian forecasting approaches.

The starting point of BPDS is the production of a prior distribution over the $m$-dimensional unknown outcome $\mathbf{y}$ – implicitly indexed by $T+1$ in a time series forecasting application – and the information set $\mathcal{H}$, encoded via the $J$ predictive models $\{h_{j}(\mathbf{y}|\mathbf{x}_{j}):1\leq j\leq J\}$, where $\mathbf{x}=(\mathbf{x}_{1},\dots,\mathbf{x}_{J})$ denotes the collection of vectors of (possibly latent) dummy variables associated with a decision. The decision maker then constructs a predictive by integrating out $\mathbf{x}$ using a ‘synthesis function’ $\alpha(\mathbf{y}|\mathbf{x})$:

The choice of the synthesis function $\alpha(\cdot|\mathbf{x})$ can be used to drive the analysis. For instance, in the case of forecast combinations, we can take $h_{j}(\mathbf{y}|\mathbf{x})=p_{j}\left(\mathbf{y}|\mathbf{x},\mathcal{M}_{j}\right)$, for some model $\mathcal{M}_{j}$, and then any set of synthesis functions $\alpha_{j}(\cdot|\mathbf{x})$ such that the combination density

is a valid density, for given weights $0<\omega_{j}<1,$ $\sum_{j=1}^{J}\omega_{j}=1.$ Specific choices of $\alpha_{j}\left(\mathbf{y},\mathbf{x}_{j}\mid\mathbf{x}\right)$ then produce different forecast combination methods (see, Johnson,, 2017, for a discussion); for example, in the case of McAlinn and West, (2019) and McAlinn et al., (2020), the synthesis function is taken to be the density of a (possibly multivariate) dynamic linear factor model.

In an attempt to ‘focus’ the BPDS approach towards decisions that are tailored to a specific user-chosen loss function underlying the analysis or decision at hand, Tallman and West, (2022) propose taking as their synthesis function, $\alpha_{j}(\mathbf{y},\mathbf{x}_{j}|\mathbf{x})=\exp\{\tau^{\prime}(\mathbf{x})S_{j}(\mathbf{y},\mathbf{x}_{j})\}$, where the score $S_{j}(\mathbf{y},\mathbf{x}_{j})$ is a $k$-dimensional vector that measures the utility one receives from realizing outcome $\mathbf{y}$ under decision $\mathbf{x}_{j}$, and $\tau(\mathbf{x})$ is a vector that weights the directional relevance of $S_{j}$$(\mathbf{y},\mathbf{x}_{j})$.

While the BPS framework, as a whole, can set the tenor of the predictions towards dynamic forecast updates that produce predictions tailored to a loss function of interest, via the choice of synthesis function $\alpha(\cdot|\cdot)$, BPS is ultimately tied to a ‘likelihood-type’ framework, or at least a log-loss function, due to the presence of the latent variables $\mathbf{x}$, which must be integrated out via assumed predictive models, $p_{j}\left(\mathbf{y}|\mathbf{x},\mathcal{M}_{j}\right)$, and with these individual predictives produced using likelihood-based Bayesian methods. While the BPDS approach can somewhat circumvent the reliance on the likelihood, due to its ability to focus on specific scores, this approach appears to be distinct from methods that entirely replace the likelihood function in the update. Therefore, a very interesting research path would involve combining the methods based on generalized posteriors discussed in Section 3.3.2 with the BPS framework.