Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
127,711 characters · 26 sections · 398 citation commands
Bayesian Forecasting in Economics and Finance: A Modern Review
\baselineskip17.5pt
The Bayesian statistical paradigm uses {the rules and language of probability} to quantify uncertainty about all unknown aspects of phenomena that generate observed data. This core characteristic of the paradigm makes it particularly suitable for forecasting, with uncertainty about the unknown values of future observations automatically expressed in terms of a probability distribution. Moreover, Bayesian methods -- in principle -- allow a user to seamlessly, and systematically, yield probabilistic forecasts that reflect uncertainty about all unknowns and that, as a consequence, condition primarily on known past events, or data: a feature that geweke2006bayesian refer to as the principle of relevant conditioning.
Indeed, the ability of Bayesian forecasters to appropriately incorporate the uncertainty associated with the production of forecasts, while utilizing all available information -- both a priori and sample information -- in a principled manner, led granger1986forecasting to conclude that:
Echoing these sentiments, in our opinion, the power of the Bayesian forecasting paradigm is a product of the paradigm's ability to treat all elements of the statistical problem necessary to produce forecasts -- future observations, past observations, parameters, latent variables, models -- as arguments of a joint probability distribution. The express probabilistic formulation of these elements, in turn, allows a Bayesian to invoke the standard rules of probability to produce a distribution for an unknown future value that is conditioned on the known past data, and is marginal of other arguments that are inherently unknown.
While this ability to marginalize all unknowns through probability calculus is the hallmark of the Bayesian approach, the benefits of the paradigm, and what ultimately in our opinion defines a `good Bayesian', is the attention to detail necessary to successfully implement Bayesian methods. In Bayesian forecasting, before we ever attempt to produce a forecast, we must first carefully enumerate all possible sources of uncertainty -- including, where possible, the set of alternative forecasting models; and construct reasonable prior beliefs for these quantities, which often include (possibly several layers of) latent variables that have a specific and delicate interaction with the observed data; always taking great care to ensure that these prior beliefs do not conflict with the observed data. Then and only then can we `turn the Bayesian crank'\ to produce the joint posterior distribution over all unknown quantities (including future values), and ultimately integrate out the quantities we are not interested in to obtain the (posterior) predictive distribution for the future values of our random variables of interest. The attention to detail necessary to produce Bayesian forecasts aims to reduce the number of implicit maintained assumptions, and what explicit assumptions are maintained (e.g. the conditioning on a particular model, or finite model set) can often be rationalized/tested against the data.
Consistent with the internal coherence of the Bayesian statistical paradigm, the basic manner in which all Bayesian forecasting problems are framed is the same. What differs however, from case to case, is the way in which the problem is solved -- i.e. the way in which the forecast distribution is accessed. To understand why this is so, it is sufficient to recognize that virtually all Bayesian quantities of interest, including forecast distributions, can be expressed as expectations of some sort. For most models that are used to predict empirically relevant data these expectations are not available in closed form. Hence, in any practical problem, implementation of Bayesian forecasting is both model- and data-dependent, and relies on advanced computational tools. Different forecasting problems -- defined by different forms and `sizes' of models and data sets -- require, in turn, different approaches to computation. The evolution of the practice of Bayesian forecasting has, as a consequence, gone hand-in-hand with developments on the computational front; with increasingly large and complex models rendered amenable to a Bayesian forecasting approach via access to modern techniques of computation.
In this review, we give a modern take on the current landscape of Bayesian forecasting. Whilst excellent textbook treatments of Bayesian forecasting are given in geweke2005contemporary and west2006bayesian, and with geweke2006bayesian reviewing specific aspects of Bayesian forecasting in a slew of practical settings, the field has advanced by leaps and bounds in the last twenty years. Therefore, we believe the time is ripe to consider a review of the subject that touches on many of the novel and exciting areas now being explored. The methodological advances we review have general applicability to all discipline areas; nevertheless, due to our own interests, expertise and experience -- and to keep the scope of the paper manageable -- we have chosen to focus primarily on applications in the economic sciences. Whilst the paper is not designed to be a treatise on Bayesian computation, sufficient details are provided to enable the practitioner to understand why numerical tools are needed in most forecasting settings, and how they are used.
The general structure of the paper is as follows. In Section (ref) we provide a short tutorial on Bayesian forecasting. This begins with an outline of the Bayesian forecasting method, followed by an overview of the computational techniques used to implement the method. In Section (ref) we then take the reader on a potted chronological tour of Bayesian forecasting, up to the present day. We begin by giving a snapshot of the forecasting problems tackled during the last decade of the 20th century (and the early years of the 21st), and the computational solutions that were adopted then -- most notably, Markov chain Monte Carlo (MCMC) algorithms. We then look at the types of `intractable' forecasting problems that are increasingly encountered in the 21st century, and provide an overview of the new computational solutions that have been proposed to tackle such problems. We also outline very recent developments in which misspecification of the forecasting model is explicitly acknowledged, and conventional likelihood-based Bayesian forecasting eschewed as a consequence; with problem-specific measures of forecast accuracy (or forecast loss) used, instead, to drive the production of forecast distributions. Section (ref) then provides the reader with more detailed reviews of contemporary Bayesian forecasting in the following four broad fields: macroeconomics, finance, marketing, and electricity pricing and demand. Section (ref) closes the paper with a brief summary of the current state of play.
Before proceeding further, we make a note about scope and language. To render the scope of the paper manageable we focus primarily on Bayesian forecasting in `time series models' -- i.e. models for random variables that are indexed by time -- and on using such models to say something about the values that these random variables will assume in the future. These future values may be informed only by past observations\ on the variable, or may also depend on the known values of covariates, or regressors. We also follow the convention in the Bayesian literature by using the terms `forecast' and `prediction' (and all of their various grammatical derivations) synonymously and interchangeably in this case, for the sake of linguistic variety. The fundamental principles of Bayesian prediction apply equally to data indexed by something other than time. The term `forecast' is not used in this case as it is a term reserved for temporal settings. The main exceptions to our focus on time series models, and forecasting per se, occur in Section (ref), in which models for cross-sectional data are used to predict customer choice in marketing settings, and Section (ref), in which models for electricity demand that have a spatial dimension are referenced.
For the sake of illustration, we assume a scalar random variable $y_{t}$, and define the $(T\times1)$ vector of observations on $y_{t}$ as $\mathbf{y} _{1:T}=(y_{1},y_{2},...,y_{T})^{\prime}$. We assume (for the moment) that $\mathbf{y}_{1:T}$ has been generated from some parametric model with likelihood $p(\mathbf{y}_{1:T}|\boldsymbol{\theta})$, with $\boldsymbol{\theta }=(\theta_{1},\theta_{2},...,\theta_{p})^{\prime}$ ${\in\Theta\subseteq}$ $\mathbb{R}^{p}$ a $p$-dimensional vector of unknown parameters, and where we possess prior beliefs on $\boldsymbol{\theta}$ specified by $p(\boldsymbol{\theta})$. Using the same symbol $\mathbf{y}_{1:T}$ to denote both the vector of observed data and the $T$-dimensional vector random variable, we define the joint distribution over $\mathbf{y}_{1:T}$ and $\boldsymbol{\theta}$ as $p(\mathbf{y}_{1:T},\boldsymbol{\theta}).$ Application of the standard rules of probability to $p(\mathbf{y} _{1:T},\boldsymbol{\theta})$ yields Bayes theorem (or Bayes rule),
where $p(\mathbf{y}_{1:T})=\int_{{\Theta}}p(\mathbf{y}_{1:T} |\boldsymbol{\theta})p(\boldsymbol{\theta})d\boldsymbol{\theta}$. Bayes theorem provides a representation for the posterior probability density function (pdf) for $\boldsymbol{\theta}$, $p(\boldsymbol{\theta }\mathbf{|y}_{1:T})$, as proportional to the product of the likelihood function and the prior. The term $p(\mathbf{y}_{1:T})$ defines the marginal likelihood, and the scale factor $\left[ p(\mathbf{y}_{1:T})\right] ^{-1}$ in ((ref)) ensures that $p(\boldsymbol{\theta}|\mathbf{y} _{1:T})$ integrates to one.
Now, define $y_{T+1}$ as the (one-step-ahead) future random variable, where we focus on one-step-ahead forecasting in Sections (ref) and (ref) merely to simplify the exposition. Assuming $y_{T+1}$ to be a continuous random variable (again, for illustration), standard probability manipulations lead to the following expression for the forecast (or predictive) pdf for $y_{T+1}:$
When no confusion arises, we also refer to $p(y_{T+1}|\mathbf{y}_{1:T})$, albeit loosely, as the forecast (or predictive) distribution, or simply as the `predictive'.\footnote{We note that $p(y_{T+1}|\mathbf{y} _{1:T})$ is sometimes referred to as a `posterior' predictive in the literature, given that it is produced by averaging the conditional predictive, $p(y_{T+1}|\boldsymbol{\theta}\mathbf{,y}_{1:T})$, with respect to the posterior density, $p(\boldsymbol{\theta}\mathbf{|y}_{1:T}).$ We do not adopt this expression, leaving it to the context to make it clear as to whether the term `predictive' is being used to refer to the distribution that is marginal of $\boldsymbol{\theta}$, $p(y_{T+1}|\mathbf{y}_{1:T})$, or that which is conditioned on $\boldsymbol{\theta}$, $p(y_{T+1}|\boldsymbol{\theta }\mathbf{,y}_{1:T})$. We also streamline the exposition by not using explicit notation for any observed covariates on which the model for $y_{t}$ may depend, and on which the predictive for $y_{T+1}$ would condition, unless this is essential.} The density $p(y_{T+1}|\mathbf{y}_{1:T})$ summarizes all uncertainty about $y_{T+1}$, conditional on the assumed model -- which underpins the structure of both the conditional predictive, $p(y_{T+1} |\boldsymbol{\theta}\mathbf{,y}_{1:T})$, and the posterior itself -- and the prior beliefs that inform $p(\boldsymbol{\theta}\mathbf{|y}_{1:T}).$ Point and interval predictions of $y_{T+1}$, and indeed any other distributional summary, can be extracted from ((ref)). In the case where the model itself is uncertain, and a finite set of parametric models, $\mathcal{M}_{1}$, $\mathcal{M}_{2}$,...,$\mathcal{M}_{K},$ is assumed to span the model space, a `model-averaged' predictive (e.g. raftery:madigan:hoeting:1997, Section 2), $p_{MA}(y_{T+1}|\mathbf{y}_{1:T})$, is produced as
where $p(y_{T+1}|\mathbf{y}_{1:T},\mathcal{M}_{k})$ denotes the density in ((ref)), but now conditioned explicitly on the $kth$ model in the set. The $kth$ posterior model probability, $p(\mathcal{M}_{k}|\mathbf{y} _{1:T})$, $k=1,2,...,K,$ is computed via a further application of Bayes theorem in which the (initial) joint distribution of interest is defined over both the model space and the space for the parameters of each of the $K$ models. Standard manipulations lead to
where
for each $k=1,2,...,K$, with $\boldsymbol{\theta}_{k}$ denoting the parameter set for the $kth$ model.
As is clear, analytical evaluation of $p(y_{T+1}|\mathbf{y}_{1:T})$ in ((ref)) requires, at the very least, a closed-form expression for $p(\boldsymbol{\theta}\mathbf{|y}_{1:T})$. Typically, however, such an expression is not available, with most posteriors being known only up to a constant of proportionality, as
The main exceptions to this occur when $p(\mathbf{y}_{1:T}|\boldsymbol{\theta })$ is from the exponential family, and either a natural conjugate, or convenient noninformative prior is adopted; specifications which may be suitable for some simple (and low-dimensional) empirical problems, but are certainly not broadly applicable in practice. Analytical evaluation of $p_{MA}(y_{T+1}|\mathbf{y}_{1:T})$ in ((ref)) also requires a closed-form expression for each $p(\mathbf{y}_{1:T}|\mathcal{M}_{k})$ (with normalization of $p(\mathcal{M}_{k}|\mathbf{y}_{1:T})$ then straightforward); once again a rare thing beyond the exponential family (and standard prior) setting. Hence the need for numerical computation to implement Bayesian forecasting in virtually all realistic empirical problems.\footnote{Numerous textbook illustrations of the material in this section can be found. In addition to geweke2005contemporary and west2006bayesian as cited above, some examples are zellner:1971, koop2003bayesian and robert:2007. We also refer the reader to steele:2020 for a recent review of Bayesian model averaging in economics.}
The form of ((ref)) makes it clear that the Bayesian predictive pdf, $p(y_{T+1}|\mathbf{y}_{1:T})$, is nothing more than the posterior expectation of the predictive that conditions on $\boldsymbol{\theta }$. Hence, accessing $p(y_{T+1}|\mathbf{y}_{1:T})$ amounts to the evaluation of an expectation. This insight is helpful, as it enables us to see many of the computational methods that are used to access $p(y_{T+1}|\mathbf{y} _{1:T})$ -- in cases where it is not available in closed form -- simply as different ways of numerically estimating an expectation.
It is convenient to group Bayesian computational methods into three categories: (i) deterministic integration (or quadrature) methods (davis1975numerical; naylor:smith:1982); (ii) exact simulation methods; and (iii) approximate methods. Given that the production of $p(y_{T+1}|\mathbf{y}_{1:T})$ involves integration over $\boldsymbol{\theta }$, only in very low-dimensional models is (i) a feasible computational approach on its own, due to the well-known `curse of dimensionality' that characterizes numerical quadrature. Hence, the computational methods in (ii) and (iii) are those most commonly adopted, and will be our focus here; noting that quadrature may play still a limited role within these alternative computational frameworks.
The methods in (ii) use simulation\ to produce $M$ draws of $\boldsymbol{\theta}$, $\boldsymbol{\theta}^{(i)}$, $i=1,2,...,M$, from the posterior $p(\boldsymbol{\theta}\mathbf{|y}_{1:T})$, which, in turn, define $M$ conditional predictives, $p(y_{T+1}|\boldsymbol{\theta}^{(i)} \mathbf{,y}_{1:T})$, $i=1,2,...,M$, the mean of which is used to estimate ((ref)). Alternatively, if it is easier to simulate from $p(y_{T+1}|\boldsymbol{\theta}^{(i)}\mathbf{,y}_{1:T})$ than to evaluate it at any point in the support of $y_{T+1}$, $M$ draws of $y_{T+1}$, $y_{T+1}^{(i)} $, $i=1,2,...,M$, are taken, one for each draw $\boldsymbol{\theta}^{(i)}$, and kernel density estimation methods used to produce an estimate of $p(y_{T+1}|\mathbf{y}_{1:T})$. Different simulation methods are distinguished by the way in which the posterior draws are produced. Methods in (ii) include Monte Carlo sampling (metropolis:ulam:1949), importance sampling (IS) (hammersley:handscomb:1964; kloek1978bayesian; geweke:1989) and MCMC sampling -- including Gibbs sampling (geman:1984; gelfand:smith90) and Metropolis-Hastings (MH) algorithms (metropolis:1953; hastings:1970) -- with MCMC being by far the most common simulation method used to compute forecast distributions in practice. The term `exact' arises from the fact that, under appropriate conditions (including convergence of the Markov chain to $p(\boldsymbol{\theta}\mathbf{|y}_{1:T})$ in the case of the MCMC algorithms), such methods all produce a $\sqrt{M}$-consistent estimate of the ordinate $p(y_{T+1}|\mathbf{y}_{1:T})$, at any point in the support of the random variable $y_{T+1}$; this estimate can thus be rendered arbitrarily accurate, for large enough $M.$
We refer the reader to: chib_greenberg_1996 and geyer2011introduction for reviews of MCMC sampling; casella:george:1992 and chibandgreenberg:1995 for descriptions of the Gibbs and MH algorithms (respectively) that are useful for practitioners; and andrieu2004computational, robert:casella:2011 and martin2022history for historical accounts of MCMC sampling. geweke2006bayesian also serves as an excellent reference on the use of these computational methods in a forecasting context. Given the critical role played by MCMC methods in the production of Bayesian forecasts, the basic principles of the algorithms are also outlined below in Section (ref); with more recent developments of both IS and MCMC -- most notably sequential Monte Carlo (SMC) (gordon:salmon:smith:1993; Chopin2020) and pseudo-marginal MCMC (beaumont:2003; andrieu:roberts:2009; andrieu:doucet:holenstein:2010)\ -- discussed briefly in Section (ref).
The methods in (iii) replace $p(\boldsymbol{\theta}\mathbf{|y} _{1:T})$ in the integrand of ((ref)) with an approximation of some sort, and evaluate the resultant integral. In so doing, such methods do not aim to estimate $p(y_{T+1}|\mathbf{y}_{1:T})$ itself, but some representation of it, defined as the expectation of $p(y_{T+1} |\boldsymbol{\theta}\mathbf{,y}_{1:T})$ with respect to the relevant posterior approximation. The methods in (iii) have been based on the principles of approximate Bayesian computation (ABC) (marin:pudlo:robert:ryder:2011; sisson2011likelihood; sisson2018handbook), Bayesian synthetic likelihood (BSL) (price2018bayesian), variational Bayes (VB) (blei2017variational), and integrated nested Laplace approximation (INLA) (rue:martino:chopin:2009), and produce what are termed `approximate' forecast, or predictive distributions. Suffice to say that the principle adopted for estimating the `approximate predictive' so defined is typically one and the same: draws of $\boldsymbol{\theta}$ from the approximate posterior (however produced) are used to produce either a sample mean of conditional predictives, or $M$ draws of $y_{T+1}$ from $p(y_{T+1} |\boldsymbol{\theta}\mathbf{,y}_{1:T})$, with kernel density estimation then applied.
The production of ((ref)) requires the computation of each model-specific predictive, plus the computation of each ((ref)). The first set of $K$ computations would proceed via the sorts of steps outlined above. Computation of the $K$ marginal likelihoods could also be performed via one of the three broad methods listed above (in particular (ii) or (iii)); however, the fact that each ((ref)) is a prior, rather than a posterior expectation does have implications for precise manner in which computation is implemented. (See ARDIA2012, and llorente2021marginal, for details).
As is clear from the brief synopsis above, it is simulation that is key to computing forecast distributions when they are not available in closed form. While the use of simulation to compute statistical quantities of interest was known by the 1970s (metropolis:ulam:1949; metropolis:1953; hammersley:handscomb:1964; hastings:1970), the technology required to perform simulation in a convenient and timely fashion was not yet available, and simulation-based computation thus remained largely out of reach. To quote geweke2006bayesian:
In the latter part of the 20th century, things changed. The increased speed and availability of desktop machines (ceruzzi:2003), allied with critical advances in simulation methodology, led to a proliferation of methods for accessing $p(y_{T+1}|\mathbf{y}_{1:T})$ via the simulation of draws from $p(\boldsymbol{\theta}\mathbf{|y}_{1:T})$. To this end, we give a brief outline of the pre-eminent posterior simulation algorithms of the 1990s (and the early 2000s): Gibbs sampling (Section (ref)), MH-within-Gibbs sampling (Section (ref)), and (MH-within-) Gibbs sampling allied with data augmentation (Section (ref)); touching on the types of forecasting models that were able to be treated via such methods, most notably the ubiquitous state space models that underpin much modern Bayesian forecasting. To keep the exposition concise, we place all algorithmic details in Appendix (ref), and reference specific algorithms from Appendix (ref) at suitable points in the text.
As a general rule, if $p(\boldsymbol{\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(\boldsymbol{\theta}\mathbf{|y}_{1:T})$ can be decomposed into recognizable densities, from which computer simulation is feasible. IS (kloek1978bayesian; geweke:1989), via use of an `importance' or `proposal' density, $q(\boldsymbol{\theta}\mathbf{|y} _{1:T})$, that matches $p(\boldsymbol{\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(\boldsymbol{\theta}\mathbf{|y}_{1:T})$ when the dimension of $\boldsymbol{\theta}$ is large, due to the difficulty of finding a $q(\boldsymbol{\theta}\mathbf{|y}_{1:T})$ that is a `good match' to $p(\boldsymbol{\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 (geyer2011introduction), draws $\boldsymbol{\theta}^{(i)}$, $i=1,2,...,M$, produced via iterative sampling from the full conditionals converge in distribution to $p(\boldsymbol{\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\textbf$y_{T+1}$\textbfin the manner described in Section (ref). Decisions about how to partition, or `block'\textbf$\boldsymbol{\theta}$\textbfneed to be made (liu:won:kon94; roberts:sahu:1997), with a view to increasing the `efficiency' of the chain which, in effect, amounts to ensuring an accurate estimate of\textbf$p(y_{T+1}|\mathbf{y}_{1:T})$ for a given number of draws,\textbf$M.$ (See Algorithm (ref) in Appendix (ref).)
CHIB1993 and McCulloch1994 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(\boldsymbol{\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:1953 \footnote{For example, dongarra2000guest rank the Metropolis algorithm proposed in metropolis:1953 as one of the 10 algorithms \textquotedblleft with the greatest influence on the development and practice of science and engineering in the 20th century\textquotedblright.}. 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.\footnote{Any non-standard probability distribution can, in principle, be drawn from using ICDF sampling. The term `Griddy Gibbs' sampling was first used by ritter: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 Bauwens1998 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(\boldsymbol{\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.\footnote{Moreover, 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 $\boldsymbol{\theta }^{(i)}$, $i=1,2,...,M$, from the MH-within-Gibbs algorithm converge in distribution to $p(\boldsymbol{\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 (ref) in Appendix (ref).)
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 (ref).
For many empirical problems in economics and related fields, a suitable model can be partitioned into two sets: static unknowns $\boldsymbol{\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( \boldsymbol{\theta},\mathbf{z} _{1:T}\right) $ is often referred to as `data augmentation', in the spirit of tanner87, 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 ((ref)) and ((ref)), 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(\boldsymbol{\theta},\mathbf{z}_{1:T} |\mathbf{y}_{1:T})$ and, thus, from $p(\boldsymbol{\theta}|\mathbf{y}_{1:T})$; an insight obtained independently by carter:kohn:1994 and fruhwirth-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(\boldsymbol{\theta}|\mathbf{z} _{1:T},\mathbf{y}_{1:T})$ and $p(\mathbf{z}_{1:T}|\boldsymbol{\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 carlin:polson:stoffer:1992, jacquier94, shepard97, kim1998svl, chib:2002, stroud2003, CHIB2006, STRICKLAND2006, ChibOMORI2007 and STRICKLAND2008. The reviews of fearnhead2011mcmc and giordani2011bayesian provide more detailed accounts and extensive referencing of this earlier literature.\footnote{We also note here the work of chib94, in which the state space representation of an autoregressive moving average (ARMA(p,q)) model (Harvey1981) was exploited, and the principle of data augmentation invoked, in order to enable an MH-within-Gibbs scheme to be applied.} (See also Appendix (ref).)
To conclude, and once again using the generic notation in ((ref)) and ((ref)), once draws have been produced from $p(\boldsymbol{\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},\boldsymbol{\theta})$ and $p(y_{T+1}|z_{T+1},\boldsymbol{\theta }\mathbf{,y}_{1:T})$, or by averaging the conditional predictives over all draws of $z_{T+1}$ and $\boldsymbol{\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 (ref). 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(\boldsymbol{\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, martin2021approximating, for a list of examples); however, particularly pertinent ones to mention here are continuous time models in finance with unknown transition densities gallant1996moments, $\alpha $-stable models for financial returns (and/or their volatility) (peters2012likelihood; {martin2019auxiliary}), and stochastic dynamic equilibrium models in economics (calvet2015accurate). 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 (tavare:balding:griffith:donnelly:1997; {rue:martino:chopin:2009}; braun2010variational; lintusaari2017fundamentals; betancourt:2018; johndrow2019mcmc). 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}|\boldsymbol{\theta})$ at each draw of $\boldsymbol{\theta}$, thereby inducing an $O(T)$ computational burden at each iteration in an MCMC chain.\footnote{We 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' bardenet2017markov. 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 (ref) do so whilst preserving the `exact' nature of the estimate of $p(y_{T+1} |\mathbf{y}_{1:T})$, whilst those in Section (ref) 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 robert2018accelerating, and martin2022history) can be summarized as follows: the use of more geometric information about the target posterior, most notably the use of Hamiltonian updates (neal:2011; hoffman2014no); the use of better MH candidate, or {proposal distributions, including those that `adapt' to previous draws (nott2005adaptive; roberts2009examples)} ; various types of combinations of multiple chains (jacob:robert:smith:2010; neal2011mcmc; neiswanger:wang:xing:2013; glynn:rhee:2014; {huber2016perfect}; jacob2019unbiased); or {the use of ex-post variance reduction methods (craiu:meng:2005; douc:robert:2011; owen2017statistically; baker2019). We refer the reader to} greenetal2015, robert2018accelerating and dunson2019hastings for detailed reviews of modern developments in MCMC, and to Jahan2020 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:salmon:smith:1993, SMC methods have evolved into a larger suite of methods used to perform both sequential and non-sequential tasks (naesseth2019elements; Chopin2020). 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:doucet:holenstein:2010; flury_shephard_2011; pitt2012some; doucet2015efficient; deligiannidis2018correlated). 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(\boldsymbol{\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 (ref) in Appendix (ref).\footnote{PMMH 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 (bardenet2017markov; quiroz2018speeding; quiroz2019speeding) 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(\boldsymbol{\theta }|\mathbf{y}_{1:T})$. Denoting the posterior approximation generically by $g(\boldsymbol{\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 ((ref)), the approximate predictive can be represented as
Given draws of $\boldsymbol{\theta}$ from $g(\boldsymbol{\theta}|\mathbf{y}_{1:T})$, and given an appropriate forward-filtering algorithm to draw from $p(z_{T}|\boldsymbol{\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\textbf$\boldsymbol{\theta}$\textbf(and\textbf$z_{T+1}$), or by applying kernel density techniques to the draws of\textbf$y_{T+1}$\textbffrom the conditional predictive.
With reference to the taxonomy of intractable problems delineated in Section (ref), the different methods of producing $g(\boldsymbol{\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)\textit{. }Both ABC and BSL avoid the need to evaluate the DGP and, hence, are feasible methods in the\textbfdoubly-intractable\textitsettings of category (a). In brief, both methods require only \textit{simulation}, not \textit{evaluation}, of the DGP. The \textit{approximation} of $p(\boldsymbol{\theta}|\mathbf{y}_{1:T})$ arises, primarily, from the fact that both methods -- in different ways -- degrade the information in the full data set,\textbf$\mathbf{y}_{1:T}$,\textbfto the information contained in a set of summary statistics, $\eta(\mathbf{y}_{1:T} )$.\textbfAs such,\textbfthe target becomes the so-called `partial' posterior for $\boldsymbol{\theta}$,\textbfwhich conditions on\textbf$\eta(\mathbf{y}_{1:T})$, rather than\textbf$\mathbf{y}_{1:T}$.\textbfThe 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 (ref) (Appendix (ref)) and (ref) (Appendix (ref)) respectively.
In contrast to ABC and BSL, VB and INLA still target the exact posterior $p(\boldsymbol{\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 ($\boldsymbol{\theta}$ in our notation). Adopting the technique of the calculus of variations, VB produces an approximation of $p(\boldsymbol{\theta}\mathbf{|y}_{1:T})$ that is `closest' to $p(\boldsymbol{\theta}\mathbf{|y}_{1:T})$ within a chosen variational family, whilst INLA applies a series of nested Laplace approximations (laplace:1774; tierney:kadane:1986; tierney:kass:kadane:1989) to a high-dimensional latent Gaussian model to produce an approximation of $p(\boldsymbol{\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(\boldsymbol{\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 (ref) and (ref) respectively.
We refer the interested reader to martin2021approximating 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. martin2021approximating 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 (ref).
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 ((ref)) -- 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: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.
loaiza2019focused{ 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 }}({gneiting2007strictly}) 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 gneiting2007strictly, 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 {gneiting2007strictly}, {Opschoor2017, and martin2022optimal} 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} (\boldsymbol{\theta}):=\sum_{t=0}^{T-1}S[p(y_{t+1}|\boldsymbol{\theta },\mathbf{y}_{1:t}),y_{t+1}]$, where $p(y_{t+1}|\boldsymbol{\theta} ,\mathbf{y}_{1:t})$ is the pdf associated with a given $P.$ Adopting the exponential updating rule proposed by bissiri:etal:2016 (see also giummole2017objective, holmes2017assigning, guedj2019, lyddon2019general, and syring2019calibrating), loaiza2019focused define the generalized (or Gibbs) posterior:
for some learning rate $\omega\geq0$, calibrated in a preliminary step. This posterior explicitly places high weight on -- or focuses on -- values of $\boldsymbol{\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}|\boldsymbol{\theta }\mathbf{,y}_{1:T})$, is misspecified, ((ref)) will -- out-of-sample -- often outperform, in the chosen rule $S(\cdot,\cdot)$, the likelihood (or {log-score})-based predictive in ((ref)), and this is demonstrated in loaiza2019focused 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 {frazierloss, in which the principles delineated here are extended to high-dimensional models, and approximations to both }$\pi_{w}(\boldsymbol{\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 loaiza2019focused and frazierloss 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:etal:2016, holmes2017assigning, lyddon2019general and syring2019calibrating.
The predictive distributions within the `plausible class' referenced above {may characterize a single dynamic structure} depending on a vector of unknown parameters, $\boldsymbol{\theta}$, {or may constitute} weighted combinations of predictives from distinct {models}, in which case $\boldsymbol{\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\footnote{See {Hall2007, ranjan2010, Geweke2011 and gneiting:2013} for early contributions to the frequentist forecast combination literature, and Wang2022 for a recent review. We note that whilst {Geweke2011 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 (Dawid1982; Dawid1985; gneiting2007probabilistic), {include }{Billio2013}, {{casarin2015jss}, casarin2015, {casarin2016}, {Pett2016}}, aastveit2018, {{{Bassetti2018, BASTURK2019}} and {CASARIN2023. Once again adopting the language of bernardo: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. Billio2013, and CASARIN2023) or implicitly allow for such a possibility (e.g. loaiza2019focused). 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 Koop2012, page 875 for an illustration of this). We refer the reader to raftery:2010 for the initial proposal of DMA, Koop2012 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.\footnote{We also refer the reader to green:1995, madigan: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 green03, marin:mengersen:robert:2005, chib2011introduction and fan2011reversible 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) (johnson2017bayesian; mcalinn2019dynamic; mcalinn2020multivariate; Aastveit2022), recently expanded to Bayesian predictive decision synthesis (BPDS) by tallman2022bayesian. 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})$: \[ p(\mathbf{y}|\mathcal{H})=\int_{\mathcal{X}}\alpha(\mathbf{y}|\mathbf{x} )\prod_{j=1}^{J}h_{j}(\mathbf{y}|\mathbf{x}_{j})d\mathbf{x}_{1}\dots d\mathbf{x}_{J}. \] 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 \[ p(\mathbf{y}|\mathcal{H})=\int_{\mathcal{X}}\frac{\sum_{j=1}^{J}\omega _{j}\alpha_{j}\left( \mathbf{y},\mathbf{x}_{j}\mid\mathbf{x}\right) p_{j}\left( \mathbf{y}\mid\mathbf{x},\mathcal{M}_{j}\right) }{\sum_{k=1} ^{J}\omega_{k}\alpha_{k}\left( \mathbf{y},\mathbf{x}_{k}\mid\mathbf{x} \right) }d\mathbf{x}_{1}\dots d\mathbf{x}_{J}. \] 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, johnson2017bayesian, for a discussion); for example, in the case of mcalinn2019dynamic and mcalinn2020multivariate, 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, tallman2022bayesian 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 (ref) with the BPS framework.
Having established the necessary details regarding the production of Bayesian forecasts in general contexts, we now review how this general probabilistic mechanism is employed to produce Bayesian forecasts in several important empirical fields. In order to produce a comprehensive and up-to-date review of each area, a range of discipline experts have been invited to write the various sections, with the authorship flagged in the section headings. This means that the style of coverage differs somewhat across sections, as suits the topic, and as fits with the perspective of the authors. However, we have aimed to retain notation that (as far as possible) is both consistent across sections and consistent with the notation used in the earlier parts of the paper and in the technical appendix, plus to ensure that the basic layout of all sections is the same. As noted earlier, other than in Section (ref) -- in which cross-sectional consumer choice data is modelled -- and in Section (ref), in which spatial models are briefly referenced, time series problems and forecasting are the primary focus.
Central banks and other policy institutions routinely collect vast amounts of time series data on key macroeconomic outcomes. One stylized fact is that these data sets often display substantial co-movements and this calls for modeling all these series jointly to produce accurate point and density forecasts. This, however, leads to large-scale models that are prone to overfitting, ultimately resulting in weak out-of-sample forecasting performance. This helps explain the popularity of Bayesian methods for macroeconomic forecasting. They can easily handle many parameters and, through appropriate prior choice, deal effectively with questions related to model and specification uncertainty in macroeconomic settings.
At a high level of generality, there are two modelling approaches used by macroeconomic forecasters. The first uses reduced-form models and imposes relatively little economic structure on the data. The second uses structural models such as dynamic stochastic general equilibrium (DSGE) models that are often estimated through Bayesian techniques; see, among many others, ALV2007, SW2007AER and DS2016. However, reduced-form approaches have proved more popular and, in this section, our focus will be on them.
As stated above, macroeconomists are typically interested in modeling the joint evolution of a set of macroeconomic quantities. To set up a general framework for understanding the types of models used for forecasting, assume that an $M$-dimensional vector $\mathbf{y}_{t}$ is related to a $K$ -dimensional vector of explanatory variables $\mathbf{x}_{t}$ through
where $g:\mathbb{R}^{K}\rightarrow\mathbb{R}^{M}$ is a function and $\bm\varepsilon_{t}$ is $\mathcal{N}(\bm0_{M},\bm\Sigma_{t})$.\footnote{Note that the Gaussianity assumption is not essential; mixtures of Gaussian distributions, for example, can be used to produce flexible error distributions if deemed necessary (see, for example, CHKMP2022, and LP2022).} This general specification nests most important reduced-form models commonly used in macroeconomics and can be used to explain the main issues that arise.
For instance, if $\mathbf{x}_{t}=(\mathbf{y}_{t-1}^{\prime},\dots ,\mathbf{y}_{t-p}^{\prime})^{\prime}$ contains $p$ lags of $\mathbf{y}_{t}$, $g(\mathbf{x}_{t})=\mathbf{A}\mathbf{x}_{t}$ is a linear function with $M\times K(=Mp)$ coefficient matrix $\mathbf{A}$, and $\bm\Sigma_{t} =\bm\Sigma$ is constant over time we have a standard vector autoregressive (VAR) model. If we set $\mathbf{x}_{t}=\bm f_{t}$ with $\bm f_{t}$ denoting a set of $Q\ll M$ latent factors and $g(\bm f_{t})=\bm\Lambda\bm f_{t}$ is linear with $\bm\Lambda$ being an $M\times Q$ matrix of factor loadings and $\bm f_{t}$ evolves according to some stochastic process (such as a VAR), we end up with a dynamic factor model stock2011dynamic. Factor augmented VARs bernanke2005measuring combine a VAR with a DFM. The dependent variables in the VAR part of the model are a subset of $\mathbf{y}_{t}$ plus a small number of factors.
Traditionally, VARs and factor models have been linear and homoskedastic. But there is a great deal of empirical evidence in most macroeconomic data sets of parameter change, both in the conditional mean and the conditional variance. This can be accommodated through particular choices for $g$ and $\bm\Sigma _{t}$. For the latter, stochastic volatility processes have proved particularly popular. For the former, various parametric forms for $g$ lead to time-varying parameter VARs (TVP-VARs) which assume that the coefficients of the VAR evolve according to a random walk. But it is also worth noting that there is an increasing literature which assumes $g$ is unknown and uses Bayesian nonparametric methods to uncover its form (see, for example, kalli2018bayesian, adrian2021multimodality, and HUBER2020).
If we set $M=1$ we obtain single-equation time series regressions, which are particularly popular in inflation forecasting (e.g. based on the Phillips curve). If we additionally set $x_{t}=1$ and allow for time-varying parameters, we can obtain models such as the unobserved components stochastic volatility (UCSV) model of SW2007 that is commonly used to forecast inflation (for recent applications, see chan2013new, stock2016core, and huber2021dynamic).
This general framework defines a class of likelihood functions. As per the outline in Section (ref), Bayesian forecasting involves multiplying a chosen likelihood function by an appropriate prior to produce a posterior which can be used to produce the predictive density. The choice of prior and computational method used for posterior and predictive inference will be case specific and we will have more to say about some interesting cases below. But a few general comments are worth noting here. First, the choice of prior matters much more in models such as the large VAR, which have a large number of parameters relative to the number observations, than in models with fewer parameters such as the UCSV model or the DFM. Second, for linear homoskedastic models with conjugate priors analytical formulae for the posterior and the one-step-ahead predictive density are available. For all other cases, MCMC methods are available. These take the general form outlined in Section (ref). However, as noted in Section (ref), MCMC methods typically do not scale well and can be computationally slow in models involving large numbers of parameters (such as large VARs) or large numbers of latent states (such as TVP-VARs). Thus, the focus of many recent papers has been on developing either improved MCMC algorithms or approximate VB methods for speeding up computation. Thirdly, our discussion so far focuses on forecasting with a single model. In practice, it is common to find that forecasts improve if many models are combined. Thus, either BMA or, alternatively, the methods outlined in Section (ref) are commonly used by macroeconomic forecasters.
With this general framework established, it is worthwhile to offer some additional detail about some of the most important 21st century developments and a discussion of how they have led to improvements in macroeconomic forecasting.
\paragraph{Large VARs}
Going back to early work such as DLS1984, Bayesian VARs have been used successfully in a variety of macroeconomic forecasting applications. Recently, they have enjoyed even greater popularity due to the rise of the large VAR. The pioneering large VAR paper was BGR2010. Subsequently, dozens of papers have used large VARs for macroeconomic forecasting (see, among many others, carriero2009forecasting, koop2013, carriero2015bayesian, GLP2015, and hauzenberger2021combining). Large VARs, involving dozens or even hundreds of dependent variables, have been found to forecast well and improve upon single-equation techniques and DFMs. Large VARs are heavily over-parameterized and, thus, Bayesian prior shrinkage has been essential in ensuring their forecasting success. We will discuss priors shortly, but at this point we highlight the fact that the use of large Bayesian VARs has been one of the major recent developments in macroeconomic forecasting.
\paragraph{Prior shrinkage in VARs}
Many different priors have been used with VARs. Traditionally, natural conjugate priors in the Minnesota tradition were used since these allowed for analytical posteriors and one-step-ahead predictives. Definitions of these priors and discussions of their properties are available in standard sources such as KK2010 and BEAR. These priors are subjective and require the user to select prior hyperparameters, most importantly those relating to the strength of prior shrinkage. In recent years, a range of alternative priors have been proposed which are more automatic, requiring fewer subjective prior choices by the researcher. For instance, GLP2015 develop methods for estimating shrinkage parameters in conjugate priors, thus avoiding the need for their subjective elicitation. Chanasy2022 also uses a conjugate prior and develops methods for selecting shrinkage parameters using a prior which relaxes some of the restrictive assumptions of the Minnesota prior. There are also a range of methods which automatically decide on the optimal degree of shrinkage for each VAR coefficient. These are the global-local shrinkage priors which are widely used with regressions and in machine learning applications, and increasingly used with VARs.\footnote{They are also used with DFMs to select the number of factors.} Global-local shrinkage priors have the form \[ a_{j}\sim\mathcal{N}(0,\psi_{j}\lambda),\quad\psi_{j}\sim f_{1},\quad \lambda\sim f_{2}, \] where $a_{j}$ is the $j^{th}$ VAR coefficient, $\lambda$ controls global shrinkage since it is common to all coefficients, and $\psi_{j}$ controls local shrinkage since it is specific to the $j^{th}$ coefficient. The densities $f_{1}$ and $f_{2}$ are mixing densities and a large range of choices of them have been proposed. One choice leads to stochastic search variable selection, used with VARs in GSN2008, koop2013 and korobilis2013, and many other references. Other choices lead to the Dirichlet-Laplace prior used with VARs by KastnerHuber2021, or the {normal-gamma} and {horseshoe} priors used in huber2019adaptive and CROSS2020; and there are many others. Since these priors are Gaussian at the first layer of the hierarchy, textbook MCMC algorithms for all the VAR parameters can be easily implemented.\footnote{In large VARs with global-local shrinkage priors, MCMC methods can nevertheless be very slow, with much faster VB methods developed in GEFANG2022.}
\paragraph{Adding stochastic volatility (SV)}
The other main development that has had a tremendous impact on applied macroeconomic forecasting in the 21st century is the development of models such as VARs that incorporate parameter change and nonlinearity. Put simply, the macroeconomic world is rarely linear and homoskedastic, and models that relax these assumptions have been found to improve macroeconomic forecasting. These improvements lie not only in point forecasts, but more importantly in density forecasts. Given the increasing interest, by central banks and academics alike, in issues such as forecast uncertainty and tail risk, the fact that these new models produce more accurate predictive densities increases their value.
A popular specification for VARs {with SV involves} factorizing the error variance-covariance matrix as $\bm\Sigma_{t}=\mathbf{A}_{0}\mathbf{H} _{t}\mathbf{A}_{0}^{\prime}$ with $\mathbf{A}_{0}$ being a lower triangular matrix with unit diagonals\footnote{$\mathbf{A}_{0}$ can also be time varying.} and $\mathbf{H}_{t}=\text{diag}(e^{h_{1t}},\dots,e^{h_{Mt}})$ being a diagonal matrix with log-volatilities evolving according to simple stochastic processes such as independent random walks or AR(1) processes. In an important contribution, clark2011real considers a VAR-SV and finds it to produce accurate point and density forecasts relative to homoskedastic models, with gains being particularly pronounced using forecast metrics involving the entire predictive density. Building on this insight, several other researchers have analyzed the role of heteroskedasticity in macroeconomic forecasting in VARs (see, for example, clark2015macroeconomic, and chiu2017forecasting) and confirm the result that using SV pays off when the focus is on obtaining accurate density forecasts. However, a problem with the standard SV specification is that the computational burden relative to homoskedastic VARs is increased enormously. This makes it difficult to do Bayesian forecasting with large VARs with SV. As a remedy, carriero2016common propose a simple common stochastic volatility (CSV) specification that assumes the shock variances to be driven by a single common volatility factor, maintaining conjugacy and thus leading to computationally efficient MCMC algorithms. They acknowledge that this model is simplistic but show that it yields much more accurate forecasts than homoskedastic VARs in a standard US macroeconomic forecasting application.
To gain more flexibility, researchers have developed algorithms that allow for estimating large VARs with $M$ independent SV processes. carriero2019large propose techniques that permit equation-by-equation estimation of such VARs and thus render estimation of larger models with SV feasible. Modified versions of this algorithm form the basis of several recent papers that combine large data sets with SV for macroeconomic forecasting (see, among others, huber2019adaptive, chan2021minnesota, and chan2021large).
\paragraph{Adding time variation in the VAR coefficients}
The previous discussion has emphasized that capturing changing error variances is key for obtaining precise forecasts. However, it may also be important to allow for structural change in the VAR coefficients themselves. One popular multivariate model that captures both changes in the VAR coefficients and error variances is the TVP-VAR-SV model proposed in primiceri2005time, which assumes that the VAR coefficients $\bm\beta_{t}=\text{vec}(\bm A_{t})$ are time-varying and evolve according to a multivariate random walk while $\bm\Sigma_{t}$ is a multivariate SV process. This model is a multivariate state space model which can be estimated using adaptations of the techniques outlined in Section (ref). The innovations to the states govern the amount of time variation in the parameters. Various shrinkage priors (often based on the global-local shrinkage priors discussed above) have been proposed that allow for a data-based decision as to whether time variation in a corresponding coefficient is necessary or not. These priors are typically elicited on the non-centered parameterization of the state space model fruhwirth2010stochastic and can help minimize overfitting concerns and produce improved forecasts.
d2013macroeconomic is an important early contribution to the macroeconomic forecasting literature using TVP models. This paper uses a small TVP-VAR with SV and shows that it produces more accurate point predictions, outperforming simpler univariate benchmarks and constant parameter VARs. One key shortcoming of this model, however, is that it only uses a small information set. This has led to several researchers proposing new methods that can be used in higher dimensions. Various approaches are possible, including models that restrict the TVP process (e.g. by imposing a factor structure, which allows for time variation in a large number of parameters to be driven by a low number of factors, see chan2020reducing). As mentioned above, shrinkage priors are used to keep the curse of dimensionality in check. These priors are typically used after transforming the model to allow for equation-by-equation estimation. Such approaches mean fairly high-dimensional TVP-VARs can be estimated without risk of over-fitting, and in a reasonable amount of time. MCMC-based forecasting with large TVP-VARs and regressions is also an active field of research and different shrinkage methods and advances in computation have led to improvements in the forecasting performance of TVP models (see, among many others, HHKO2021, and HKO2021). However, it is worth noting that if computation does become a concern, approximate methods (e.g. using the VB methods outlined in Section (ref)) can be used. Approaches which avoid the need for MCMC are developed in koop2013large and koop2018variational. In the former paper the authors propose large approximate TVP-VARs based on forgetting factors, whereas in the latter they use VB techniques to forecast inflation with large TVP regression models.
\paragraph{Bayesian nonparametric VARs}
Up to this point we have assumed that the conditional mean function $g$ takes a known form. However, it could be that the functional form is unknown. Bayesian nonparametric techniques, such as Bayesian additive regression trees chipman2010bart, Gaussian processes and kernel regressions adrian2021multimodality or infinite mixtures kalli2018bayesian, allow the researcher to uncover such unknown functional forms and produce precise macroeconomic forecasts. In general, they have had great success, but they have been found to be particularly useful in studies that focus on the tails of predictive distributions or on the handling of outliers such as the ones experienced during the pandemic (see, for example, HUBER2020, and clark2022forecasting).
kalli2018bayesian propose a nonparametric VAR that builds on an infinite mixture model with the mixture weights being driven by the lagged endogenous variables. They show, using US and UK data, that their model yields competitive forecasts, with accuracy gains in terms of point and density predictions increasing sharply for higher forecast horizons. CHKMP2022 use BART-based VARs to perform tail forecasting of US output, unemployment and inflation in {real time}, finding that nonparametric techniques work well in the tails and for {higher-order} forecasts. With a particular focus on predictive accuracy during the pandemic, HUBER2020 develop mixed frequency nonparametric VARs and show that these models yield substantially more precise nowcasts during the Covid-19 period.
\paragraph{Conclusions and further directions}
We have outlined how Bayesian methods have been used successfully for macroeconomic forecasting. Most of the discussion has related to VARs, which are a class of models where Bayesian methods have proved particularly popular. But it is worth noting that empirically-relevant extensions (e.g. SV or TVP) can be added to other multivariate time series models such as DFMs or FAVARs, as can the VAR prior shrinkage methods (e.g. global-local shrinkage methods) we have discussed. It is also worth noting that we have focused on models that do not restrict the coefficients. However, restricted VARs are often used for forecasting. For instance, vector error correction models (which impose cointegrating restrictions) or multi-country VARs such as global VARs are restricted VARs.
We have also focused on forecasting as opposed to the closely related field of nowcasting. Mixed frequency VARs, which jointly model quickly-released, high-frequency variables (e.g. monthly variables such as surveys, employment and inflation) and slowly-released, low-frequency variables (e.g. quarterly variables such as GDP), have proved very popular with nowcasters. Bayesian methods are typically used with such models (see, for example, SS2015, HUBER2020, KMMP, and stackedVAR) and, in real-time nowcasting exercises they tend to perform well.
A pertinent question in financial analysis is whether the risks associated with financial assets -- and the prices of those risks -- are predictable in ways that are useful in applications such as portfolio allocation, risk management and derivative pricing. With risk factors typically being represented as latent distributional features of observable financial variables, it follows that two key goals in the statistical analysis of financial problems are: i) The accurate prediction of latent distributional features; and ii) The development of complex, non-linear state space models to underpin this prediction.
Both of these goals lend themselves naturally to a Bayesian treatment given, in turn, the automatic production of predictive distributions via the Bayesian paradigm, and the swathe of computational methods available to estimate complex models -- most notably those with a latent variable structure. In particular, the growth in financial derivatives markets from the 1990s onwards has generated the need to model the underlying asset as a continuous time process, almost always augmented with a continuous time process for the asset volatility, and often via a jump diffusion. Such models -- whilst `convenient' in the sense of allowing for closed-form solutions for derivative prices -- are challenging from a statistical point of view, given that they typically need to be treated as a (discretized) non-linear state-space model, and may require multiple sources of data to enable separate identification of model parameters and risk premia. Estimation of and forecasting with such models is nevertheless computationally feasible via Bayesian methods, with MCMC algorithms of one form or another forming the backbone of the early treatments (eraker2001; eraker2003; eraker2004; forbes2007; Johannes2009).
We refer the reader to jacquier2011bayesian and JOHANNES2010 for comprehensive reviews of the application of Bayesian methods in finance up to the first decade of the 21st century. The coverage includes, in short, Bayesian approaches to: portfolio allocation, return predictability, asset pricing, volatility, covariance, `beta' and `value at risk' prediction, continuous time models (and discretized versions thereof), interest rate modelling, and derivative (e.g. option) pricing. Our goal in the current review is to outline the more recent advances that have evolved over the last decade, in particular those that have exploited (in one way or another) new methodological advances, new sources of data, and modern computational techniques. In order, we shall briefly review: the use of diverse data sets, including derivative prices and high-frequency measures of financial quantities; the treatment of DGPs that are unavailable in closed form; the analysis of high-dimensional models; and the application of non-parametric modelling.
\paragraph{Multiple sources of financial data}
It is now a well-established fact that the constant volatility feature of a geometric diffusion process for a financial asset price is inconsistent with both the observed dynamics in return volatility and the excess kurtosis and skewness that characterizes the typical empirical return distribution; see bollerslev:chou:kroner:1992 for an early review. The option pricing literature supports this finding, with certain empirical regularities, such as `implied volatility smiles', seen as evidence that asset prices deviate from the geometric Brownian motion assumption that underlies the black1973 option price bakshi1997,hafner2001,lim2005. Hence, the 21st century has seen the proliferation of many alternative specifications for asset prices, and associated theoretical derivative prices, most of which are nested in a general framework of (discretized) bivariate jump diffusion models for the asset itself and its volatility. Allied with these developments has been the growth in access to transaction-level `high-frequency' data -- in both the spot and options markets -- which, in itself, has spawned new approaches to inference and forecasting in the financial sphere.
The Bayesian literature has brought to bear on this problem the power of computational methods -- both established, and more recent -- to enable the multivariate state space models that have emerged from this literature to be estimated, and probabilistic predictions of all dynamic variables -- the return itself, volatility, random jumps (in either the return or the volatility, or both), and various risk premia -- to be produced. With reference to the generic notation for a state space model in ((ref)) and ((ref)), Bayesian approaches over the last decade can be categorized according to the specification adopted for the (multivariate) measurement at time $t$, $\mathbf{y}_{t}$ and, hence, for the (multivariate) state, $\mathbf{z}_{t}$, being modelled and forecast. Some work exploits data from both the spot and options market to predict volatility and its risk premia (maneesoonthorn2012), and option prices (Yu2011; carverhill2022\footnote{We note that whilst a time series model is constructed in the case of these two references, the (out-of-sample) prediction of option prices is across the cross section of strike prices and maturities. We also make note of FULOP2019 who exploit spot and options data to produce filtered estimates (as opposed to strictly out-of-sample predictions) of latent volatility and price jump intensity.}); other work combines `low-frequency' daily observations on returns with high-frequency measures of volatility and/or price jumps to predict (in some combination) returns, volatility, and the size and occurrence of price jumps jin2013modeling,maneesoonthorn2017,frazier2019approximate; whilst further work combines daily returns with futures prices in predicting various financial quantities of interest FILECCIA2018,gonzato2021.
\paragraph{Financial models that are `unavailable'}
All but one of the papers cited in the previous paragraphs share a common feature - namely, a DGP that can be expressed as a probability density (or mass) function. With reference to ((ref)), it is the availability of a closed form for $p(\mathbf{y}_{1:T},\mathbf{z}_{1:T}|\boldsymbol{\theta })=p(\mathbf{y}_{1:T}|\mathbf{z}_{1:T},\boldsymbol{\theta})p(\mathbf{z} _{1:T}|\boldsymbol{\theta})$, that renders feasible the MCMC methods used in the said works. In contrast, frazier2019approximate adopt a process for the latent log-volatility that is driven by an $\alpha$-stable innovation, such that $p(\mathbf{z}_{1:T}|\boldsymbol{\theta})$ is unavailable, and MCMC infeasible as a consequence. Instead, ABC is adopted for inference, and an approximate predictive of the form of ((ref)) produced instead. In addition to providing theoretical validation of the approach, the authors demonstrate, in range of different simulation settings, that despite inaccuracy at the posterior level, the approximate predictive is always a very close match to the exact predictive. Related work in which an ABC method is used to conduct forecasting appears in canale2016, konkamking2019, VIRBICKAITE2020 and pesonen2022. ABC treatment of a conditional likelihood for a time series of financial returns, $p(\mathbf{y}_{1:T}|\mathbf{z}_{1:T},\boldsymbol{\theta})$, that is unavailable in closed form is also investigated in CREEL2015, martin2019auxiliary and chakraborty2022modularized, with chakraborty2022modularized proposing a modularized version of ABC. For other recent Bayesian treatments of intractable models of this sort that continue to exploit MCMC principles (with or without an ABC component), see vankov2019filtering and muller2021estimation.\footnote{The citation of CREEL2015, martin2019auxiliary, vankov2019filtering and muller2021estimation is relevant to this review, despite these references not having an explicit component on forecasting.}
\paragraph{Large financial models}
Thus far, we have reviewed Bayesian treatments of models for single financial assets. That is, the models may have specified multiple latent components, and potentially multiple measurements, but they still aim to explain (and forecast) quantities related to a single asset. Models for multiple assets are also critically important in financial applications, with the relationship between financial assets determining the extent to which diversification can be achieved, as well as how risks permeate the various sectors of the financial market. Indeed, Bayesian methods are particularly suitable for dealing with such multivariate models, since the dimensionality of $\mathbf{z}_{1:T}$ is typically much larger than that of $\mathbf{y}_{1:T}$ and, hence, challenging to deal with via any other means.
chib2009Review provide an early review of the Bayesian analysis of multivariate SV models, with all work up to this point utilizing traditional MCMC techniques, and the statistical and predictive analysis limited to relatively low-dimensional systems (up to ten assets). Subsequent work has focused on the development of more flexible multivariate distributions nakajima2017bayesian, and the use of sparse factor structures and shrinkage priors in constructing larger-dimensional models zhou2014bayesian,kastner2017efficient,BASTURK2019. More recently, with the advances made in VB methods, inference and prediction in very large-dimensional financial models is now possible (gunawan2021variational; Chansv2022; frazierssm2021; quiroz2018gaussian; zhangw2023). There is also a growing interest in the prediction of co-movements of various sorts, with: bernardi2015bayesian predicting the interdependence between U.S. stocks with Bayesian time-varying quantile regressions; geraci2018measuring capturing and predicting the interconnectedness of financial institutions through Bayesian time-varying VARs; and alexopoulos2022bayesian modelling and predicting common jump factors in a large panel of financial returns.
\paragraph{Bayesian nonparametric modelling in finance }
As noted, simple parametric assumptions such as additive Gaussian innovations are inconsistent with the stylized features of financial data. Whilst more suitable non-Gaussian/non-linear models can be built (as highlighted above), Bayesian nonparametric modelling allows for further flexibility via the incorporation of Dirichlet process mixture (DPM) structures. Such an approach has been shown to provide robustness to distributional assumptions and can improve point forecasts, but the main gain has been significant improvements in the accuracy of predictive densities, and of risk measures derived from those densities. The advancement of the literature in this direction has been aided by the stick-breaking representation sethuraman1994constructive and the introduction of the slice sampler walker2007sampling,kalli2011slice.
jensen2010bayesian introduce an extension to a standard SV model to capture the unknown return innovation distribution via a DPM. The DPM specification has also been inserted into other popular models in finance, with: jensen2014estimating adopting a DPM to jointly model the return and future log-volatility distribution; delatola2013bayesian capturing the so-called leverage effect; AUSIN2014350 applying a DPM to univariate GARCH models; and Kalli-Griffin:2015 using Bayesian nonparametric modelling to aggregate autoregressive processes to produce an SV model with long-range dependence. Extensions to multivariate financial models have also occurred: in a multivariate GARCH setting in Jensen-Maheu:2013; and in a Cholesky-type multivariate SV model in ZTW:2020.
A potential drawback of the DPM model is that it neglects time dependence in the unknown distribution. An important extension of the DPM prior is the hierarchical Dirichlet process of teh2006hierarchical, which allows for the construction of a prior for an infinite hidden Markov model (IHMM), which allows for time dependence in a flexible manner. The introduction of the beam sampler of van2008beam, which extends the slice sampler, renders conventional posterior sampling methods for finite-state Markov switching models CHIB1996 feasible in the IHMM. The IHMM structure has been used to model the univariate GARCH distribution dufays2016infinite, and the multivariate GARCH distribution Robin:2022; and to provide a nonparametric model for realized measures, including realized covariance matrices jin2016bayesian,liu-maheu:2018,jin2019bayesian, with all papers documenting very large improvements in density forecast accuracy from the IHMM. Other applications of the IHMM include: shi-yong:2016, who use the IHMM to date and forecast speculative bubbles, and who also adopt a version with GARCH effects; yang2019, who studies the relationship between stock returns and real growth with a multivariate IHMM model; and, more recently, JIN2022302, who employ the DPM prior in the infinite Markov pooling of predictive distributions, with forecasting applications to interest rates, realized covariances and asset returns. Other approaches to time dependence in Bayesian nonparametrics for finance include griffin2011stick, who introduce a time-dependent stick breaking process in a general setting and develop an SV model for returns. More recently, SKL:2020 use a weighted DPM to forecast return distributions, while SHAMSI:2021 allows for lagged covariates to impact the weights in the DPM model through a probit stick-breaking process.
Bayesian methods are applied to a wide range of marketing problems; see rossi2003bayesian for a review of the early literature. More recently, these methods have been increasingly used for the purpose of prediction, for instance in customer choice behaviour toubia2019extracting,araya2022identifying, customer demand posch2022bayesian, customer satisfaction mittal2021improving, dynamic pricing bastani2022meta, advertising effectiveness danaher2020advertising,loaiza2022fast and recommender systems ansari2018probabilistic. Given the large variety of marketing applications, we focus in this section on the modelling of customer choice to illustrate the key principles of Bayesian prediction in marketing problems.
A common problem in marketing is that of setting the price level of a set of products so that total profits are maximized. To estimate these optimal prices, predictions of how customers will react to price changes are crucial. Predictions of customer choices under different marketing environments can be constructed by choice models. These models are estimated using data about the product choices of customers in the marketplace, a survey, an experiment, etc. rossi2012bayesian.
An example of a prediction of interest in this context is the predicted purchase probability of a customer for a particular product as a function of its own price or the price of another product. The predicted purchase probability can be constructed for a customer for which only a few choices are observed, or for a new customer for which we do not observe choices in the data.\footnote{Although this section, as noted in the Introduction, focuses on prediction using cross-sectional data, choice models can also be applied to the forecasting of future choice probabilities by using time series data mccormick2012dynamic or panel data gilbride2004choice,terui2011effect.}
The two most popular models used to predict choice behaviour are the multinomial logit and multinomial probit models. The multinomial logit model imposes the independence of irrelevant alternatives (IIA) property mcfadden1989method, which means that it cannot capture general substitution patterns among choice alternatives. The IIA property of this model can be relaxed under certain assumptions by extending the multinomial logit model to a nested logit model poirier1996bayesian,lahiri2002bayesian or\ a random parameter logit model train2009discrete.
On the other hand, the multinomial probit model does not impose the IIA property, and as such is commonly used in the analysis of economic choice behaviour, where complementary and substitution effects are important. For instance, the multinomial probit model has been recently used in the analysis of car choices karmakar2021understanding, grocery brand choices miyazaki2021dynamic, employment choices mishkin2021gender, and car parking choices paleti2018generalized. The remainder of this section presents a review of Bayesian prediction based on the multinomial probit model.
\paragraph{Multinomial probit model specification}
The variable of interest is $y_{i}\in\{0,1,2,\dots,J\}$, which indicates the choice made by individual $i$ among a set of $J+1$ alternatives. This choice is modeled to be conditional on a set of $J$ latent utilities $\mathbf{z} _{i}=\left( z_{i1},\dots,z_{iJ}\right) ^{\prime}$, so that the conditional pmf is defined as
where $p(y_{i}|\mathbf{z}_{i})=\text{Pr}(Y_{i}=y_{i}|\mathbf{z}_{i})$, $z_{iy_{i}}$ is the $y_{i}$-th element of $\mathbf{z}_{i}$, with $y_{i}>0$, and $I[A]$ is one if statement $A$ is true and zero otherwise. The base category $j=0$ is one of the choice alternatives, which is selected a priori. The base category is observed whenever all the latent utilities are less than zero.
The utilities are expressed in terms of $r$ predictors via a linear Gaussian model,
where $\phi_{J}\left( \mathbf{z};\boldsymbol{\mu},C\right) $ denotes a $J$-variate normal density with mean $\boldsymbol{\mu}$ and covariance matrix $C$, $X_{i}$ a $J\times r$ matrix of predictor values, $\boldsymbol{\beta}$ an $r$-dimensional vector of coefficients, and $\Sigma$ a covariance matrix that captures complementary and substitution effects between the choice alternatives.
Combined, (ref) and (ref) give rise to the augmented likelihood function of the multinomial probit model
where $\boldsymbol{\theta}=\{\boldsymbol{\beta},\Sigma\}$, $\mathbf{y} =\{y_{i}\}_{i=1}^{n}$, $\mathbf{z}=\{\mathbf{z}_{i}\}_{i=1}^{n}$, and $X=\{X_{i}\}_{i=1}^{n}$, with $n$ the total number of individuals. For a given prior distribution $p(\boldsymbol{\theta})$, the augmented posterior distribution of the model is given as
albert1993bayesian were the first to propose the use of data augmentation (see Section (ref) herein) for conducting Bayesian analysis of the multinomial probit model.
\paragraph{The predictive distribution}
Consider now an individual $s$, with predictor values $X_{s}$, whose choice behaviour we would like to predict. The predictive for individual $s$, can be written as
from which the predictive choice probabilities $\text{Pr}(Y_{s}=j|X_{s} ,\mathbf{y},X)=p(j|X_{s},\mathbf{y},X)$ can be constructed. The specification and computation of the predictive distribution in (ref) poses three key challenges.
First, $p({y}_{s}|\mathbf{z}_{s})$ requires a choice of base category. This choice affects the prior predictive choice probabilities, and hence the (posterior) predictive choice probabilities can be sensitive to the choice of base category; see burgette2012trace. burgette2021symmetric propose a symmetric prior specification to address this problem. The parameters $\boldsymbol{\theta}$ are not identified under this prior, but this does not affect the predicted probabilities.
Second, the parameters $\boldsymbol{\theta}$ lack scale identification, as $p(y_{i}|\mathbf{z}_{i})=p(y_{i}|c\mathbf{z}_{i})$ for any positive scalar $c$. Different solutions have been proposed to fix the scale, all based on a constraint on the specification of $\Sigma$. For instance, mcculloch2000bayesian fix the first leading element of $\Sigma$ to unity. This approach is sensitive to the ordering of the choice categories in the model. burgette2012trace fix the trace of $\Sigma$, which is invariant to the way in which the choice categories enter the model.
Third, the computation of $p({y}_{s}|X_{s},\mathbf{y},X)$ involves the evaluation of the integrals over the latent utilities in $\mathbf{z}_{s}$ and $\mathbf{z}$. Since no analytical solution for these integrals is available, they are solved with MCMC sampling steps. The latent utility of each choice category is sampled from a univariate truncated normal, conditional on the latent utilities for all the other choice alternatives, for each individual mcculloch1994exact. Conditional on the draws for the latent utilities, sampling $\boldsymbol{\beta}$ from its full conditional is straightforward. Generating from the conditional distribution of $\Sigma$ is nonstandard as the scale restrictions on $\Sigma$ have to be taken into account.
\paragraph{Scalable Bayesian prediction}
In addition to the challenges delineated above, it is difficult to scale $p({y}_{s}|X_{s},\mathbf{y},X)$ to problems with large choice sets or a large number of observations. Recent advances in the computation of the predictive have focused on tackling the scalability issues in $J$ and $n$, as we discuss below.
When considering a full covariance matrix specification for $\Sigma$, the total number of parameters increases quadratically with $J$. For problems with large choice sets and small samples, this implies that the ratio of total number of parameters to total number of observations is large, making it difficult to construct accurate predictions. loaiza2021scalable propose a spherical transformation of the covariance matrix of the latent utilities that imposes a parsimonious factor structure and a trace restriction. As a result, the total number of parameters grows only linearly with $J$. The authors demonstrate that this parsimonious structure leads to improved predictive performance over full covariance matrix specifications.
Additionally, as noted above, the construction of the predictive entails evaluation of the integral over the latent utilities $\boldsymbol{z}$. Although MCMC is able to solve this integral, it does so by generating the utility vector for each individual from a multivariate truncated normal, which is a computationally costly exercise mcculloch1994exact,botev2017normal. This renders MCMC algorithms impractical for problems where a large $n$ is considered.
VB can be employed to tackle problems with large $n$. Adapting the generic descriptions of VB\ in Section (ref) and Appendix (ref), the application of VB in this setting considers the class of approximating densities $\mathcal{Q}$ with elements $q_{{\lambda}}(\boldsymbol{\theta },\mathbf{z})\in\mathcal{Q}$, indexed by a variational parameter vector $\boldsymbol{\lambda}$. The exact augmented posterior is approximated by $q_{\hat{\lambda}}(\boldsymbol{\theta},\mathbf{z})$ with an optimal variational parameter vector equal to
where KL denotes the Kullback-Leibler divergence. The variational predictive is then constructed as \[ \hat{p}_{\lambda}({y}_{s}|X_{s},\mathbf{y},X)=\int_{\Theta}\int_{\mathbf{z} _{s}}p({y}_{s}|\mathbf{z}_{s})p(\mathbf{z}_{s}|\boldsymbol{\theta} ,X_{s})d\mathbf{z}_{s}\int_{\mathbf{z}}q_{\hat{\lambda}}(\boldsymbol{\theta },\mathbf{z})d\mathbf{z}\,d\boldsymbol{\theta}. \] Calibration of the variational approximation requires a scale-identified expression for $p(\mathbf{y},\mathbf{z}|\boldsymbol{\theta},X)$. To achieve this, girolami2006variational consider an identity matrix covariance structure, while fasano2022class fix $\Sigma$ at predetermined values. loaiza2022mnpvb propose a method for a multinomial probit model with a factor covariance structure. This method uses the hybrid variational approximation $q_{{\lambda}}(\boldsymbol{\theta},\mathbf{z})=q_{{\lambda} }(\boldsymbol{\theta})p(\mathbf{z}|\mathbf{y},\boldsymbol{\theta},X)$ introduced by loaiza2022fast.
Forecasting in electricity markets is critical for efficient day-to-day operation of power grids, long-term planning of infrastructure and increasingly, at a disaggregated level, for the management of smart grids. This section will cover forecasting electricity prices, electricity load/demand and generation by source of power, primarily wind and solar. Hereafter these problems will collectively be referred to as `electricity forecasting'. Motivations for electricity forecasting can be found in general reviews such as Wer2014 for price forecasting, LinEtal2019 for load forecasting, Ant2016 for solar power forecasting and GieKar2017 for wind power forecasting. These reviews indicate that the majority of work in electricity forecasting does not employ a Bayesian approach; however notwithstanding this, Bayesian methods have found success in the field.
There are very few instances of Bayesian forecasting in electricity markets that predate the early 2000s, although we now cover some notable exceptions. Bun1980 consider the case of updating load forecasts in an online fashion by computing a Bayesian model average of load profiles of a cloudy and a sunny day. Meanwhile, Bayesian VARs have been used by Gun1987, BecSol1994 and JouEtal1995 to forecast energy demand, nuclear power generation and demand prices and consumption respectively. A Bayesian VAR shrinks autoregressive coefficients to either a random walk or white noise depending on whether data are stationary or non-stationary and was popularized in macroeconomics by DLS1984 (see also Section (ref)). The performance of Bayesian VARs in early electricity forecasting applications is mixed; BecSol1994 find evidence in favour of Bayesian autoregression, JouEtal1995 find that Bayesian VARs are effective for forecasting demand, but not price, while Gun1987 does not find any improvement at all from using Bayesian VARs rather than conventional autoregressive integrated moving average (ARIMA) models.
With the advent and popularization of MCMC methods, Bayesian forecasting has begun to find greater success in the field of electricity forecasting. In the literature of roughly the past two decades, there are three common major motivations for using Bayesian forecasting, two of which have antecedents in the earlier literature. The first is the use of `Bayesian models'\footnote{By a `Bayesian model' we generally mean a model with a prior and likelihood estimated by Bayesian inference. Bayesian methods for finding tuning parameters such as the automatic relevance determination HipTay2010, and Bayesian optimisation lie beyond the scope of this section.}, which have now grown well beyond Bayesian VARs to include models with latent volatilities, models with a spatial dimension, and Bayesian neural networks. The second is the use of BMA for forecast combination. The third is the production of full probabilistic forecasts via Bayesian computation. These are now each discussed in turn.
\paragraph{Bayesian models}
The structure inherent in many electricity forecasting problems provides a motivation for the innovative use of priors to improve forecasting accuracy. Although the early literature cited before found somewhat ambiguous results when comparing Bayesian VARs to classical alternatives, more recent work finds evidence in favour of a Bayesian approach; see RavEtal2015 for point forecasts and GiaEtal2020 for both point and density forecasts. An important aspect of this work is the exploitation of the intraday nature of the data, since typically hourly prices are stacked in a VAR model. The intraday structure lends itself to priors that shrink parameters corresponding to consecutive hours of the day that are close to one another. An early application of this approach can be seen in CotSmi2003.
Since electricity data are increasingly available not only at a high temporal frequency but also at a high spatial resolution, there are further examples in the literature of using priors to exploit neighbourhood structure. Examples include OhtEtal2010 who use spatial ARMA processes to predict electricity load in nine Japanese regions, and GilEtal2019 who use spatio-temporal Gaussian processes to forecast residential-level electricity demand. Even where spatial information is unavailable, hierarchical models estimated using Bayesian methods have been used to produce disaggregate energy demand forecasts; examples can be found in MorNak2014 and WanEtal2017 who use Gaussian processes, and GriEtal2021 who use regression. Informative hierarchical priors have been used in instances where data sets are small in size, or unavailable; for example, PezEtal2006 elicit priors for future trajectories of temperature in the winter using past observations, and LauEtAl2015 elicit priors for the electricity demand of `non-metered' households using data on `metered' households.
While the aforementioned examples take a Bayesian approach to exploit the use of priors in novel ways, another strain of the Bayesian forecasting literature is based on estimating models with latent variables. Examples in electricity forecasting include a latent jump process for price spikes ChaEtal2014 and SV models (Smi2010; KosKos2019). Also, in recent years, Bayesian analysis of machine learning models has become increasingly popular. This includes neural network models BruEtal2019,GhaEtal2019,CapEtal2020, where VB is typically used. Also, Bayesian regression trees (see Section (ref)) have been applied to electricity forecasting by NatEtal2011 and AliEtal2019, who find that they outperform non-Bayesian counterparts. Finally, there is an extensive literature on using Bayesian networks for forecasting in energy; see AdeEtal2020 for a review of these methods in forecasting wind generation.
\paragraph{Bayesian model averaging (BMA)}
As noted earlier, the importance of forecast combination is widely appreciated in the forecasting literature. Whilst, as highlighted in Section (ref), many different Bayesian approaches to forecast combination have now been explored, BMA remains a very important method in the sphere of electricity forecasting. As described in Section (ref), BMA uses posterior model probabilities as combination weights. Whenever the choice of model is parameterized, the predictive density has an interpretation as a forecast combination. Examples include Smi2000 who combines forecasts from regression models that include different predictor sets, and PanSmi2008 who average over models with different combinations of skew and symmetric marginal distributions.
It is also common in the electricity forecasting literature to produce point forecasts from different models and then combine these using BMA as a post-processing step. This approach grew out of research combining ensembles of forecasts from numerical weather predictions (NWPs) RafEtal2005,SloEtal2010. In the NWP setting, forecasts are the outputs of deterministic physical models. Statistical models are then formed by assuming that for $k=1,\dots,K$, $p(y_{t}|a_{k},b_{k},f_{k},\sigma ^{2},\mathcal{M}_{k})\sim N(a_{k}+b_{k}f_{k},\sigma^{2})$, where $f_{k}$ is the $k^{th}$ NWP and $a_{k}$, $b_{k}$ and $\sigma^{2}$ are additional parameters. These statistical models are then combined using the usual BMA machinery described by ((ref)), with the key distinction being that posterior model probabilities are replaced with $p(M_{k}|y_{T-L+1:T})$, where $L$ is the length of the window. Uncertainty over $a_{k}$, $b_{k}$ and $\sigma^{2}$ is integrated out in the usual way, and there are no additional parameters since the $f_{k}$ are obtained deterministically. This approach has been used in energy forecasting by Coe2006, who motivate forecasting rainfall as a input into forecasting generation from hydroelectric dams, and Du2018 who uses wind forecasts to predict generation from wind farms.
The work of RafEtal2005 has been subsequently extended to the case where the forecasts $f_{k}$ are not the outputs of deterministic physical models but are point forecasts from statistical models, each with their own unknown parameters. For example NowEtal2014 adopt the approach of RafEtal2005 but where the $f_{k}$ are obtained from statistical time series models with parameters estimated using frequentist techniques. {This approach is not fully Bayesian} (despite being referred to as BMA in the literature),{ since although the model average integrates over the uncertainty in $a_{k}$, $b_{k}$ and $\sigma^{2}$ it does not integrate over uncertainty in the parameters of the underlying time series models used to generate the point forecasts $f_{k}$.\footnote{The same point does not apply when combining ensembles from NWPs since the forecasting models are deterministic.}} In a similar vein, HasEtal2015 and RazEtal2017 combine electricity load forecasts from different neural networks.
\paragraph{Probabilistic forecasting}
A common motivation for taking a Bayesian approach is the ease with which the computational machinery of MCMC or approximate methods produces a full predictive density rather than only point forecasts. Key operational decisions in electricity forecasting depend on quantities other than the predicted mean; see NowWer2018 and references therein for discussion. While the importance of probabilistic forecasting is often highlighted in Bayesian papers it is not always the case that forecasts are evaluated in a way that assesses the quality of the full predictive distribution\footnote{We note that in some cases this is challenging; for example for long-run forecasts as in daSEtal2019.}. For example, often probabilistic forecasts are summarized by prediction intervals, and the empirical coverage of these intervals used as a means of checking model quality; for an early example see PezEtal2006, and more recently WanEtal2017 and KosKos2019 , where the latter show that Bayesian methods compare favourably to non-Bayesian alternatives for forecasting electricity prices. KosKos2019 also evaluate $\alpha$-level quantile forecasts $\hat{q} _{t}$ using the pinball loss, \[ L_{\alpha}(y_{t},\hat{q}_{t})=\alpha(y_{t}-\hat{q}_{t})I[y_{t}\geq\hat{q} _{t}]+(1-\alpha)(\hat{q}_{t}-y_{t})I[y_{t}<\hat{q}_{t}]\,. \] YanEtal2019 and SunEtal2019 also use pinball loss to evaluate forecasts of residential-level load (net of solar PV generation in the latter case).
However, the use of scoring rules gneiting2007strictly and, hence, the explicit recognition of the distributional form of the forecasts, is becoming increasingly popular as a means of evaluating predictive distributions in both Bayesian and non-Bayesian electricity forecasting. The {continuously ranked probability score }(gneiting2007strictly) is particularly amenable to Bayesian inference since it is usually approximated using a Monte Carlo sample from the predictive density. For an early example of its use in Bayesian electricity forecasting see PanSmi2008; for later examples, see BraDeF2015, BruEtal2019 and GiaEtal2020. Other scoring rules are less commonly used in the Bayesian electricity forecasting literature, although OhtEtal2010, where the log score is used, is a notable exception.
Bayesian forecasting is underpinned by a single core principle: uncertainty about the future value of a random variable is expressed using a probability distribution, where the form of that distribution reflects -- in turn -- uncertainty about all other unknowns on which the investigator chooses not to condition.
While this principled approach to forecasting is arguably one of the most compelling features of the paradigm, the challenge has, potentially, been in the implementation of Bayesian forecasting: namely, computing the expectation that defines the predictive distribution, most particularly when accessing (draws from) the posterior itself is difficult. And as models have become larger and more challenging, and as data sets have grown `bigger', this problem of accessing the exact posterior has only increased. However, as this review has demonstrated, the expansion of the forecasting problems being tackled has gone hand-in-hand with the development of new and improved computational methods designed expressly to access challenging posteriors, and in a reasonable computing time. Notably, when it comes to accurate forecasting, somewhat crude approximations of the posterior have been found to still yield accurate predictions; meaning that Bayesian forecasting remains viable for large and complex models for which approximate computation of posteriors is the only feasible approach.
The more fundamental problem of model misspecification can also be managed, by moving away from the conventional likelihood-based Bayesian updating and allowing forecast accuracy itself -- and its link to the future decisions that depend on that accuracy -- to drive the updating. This, in turn, ensures that forecasts are `fit for purpose', despite the inevitable misspecification of the forecasting model. Allied with the computational power that now drives the Bayesian engine, this ability to generalize the paradigm beyond its traditional links with the likelihood principle is a potent, if not yet fully realized, force in forecasting.