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.
86,880 characters · 18 sections · 83 citation commands
Loss-Based Variational Bayes Prediction
\editor
The conventional paradigm for Bayesian prediction is underpinned by the assumption that the true data generating process is either equivalent to the predictive model adopted, or spanned by a finite set of models over which we average. Of late however, recognition of the unrealistic nature of such an assumption, allied with an increased interest in driving prediction by problem-specific measures of accuracy, {or loss, }have led to alternative approaches. Whilst antecedents of these new principles are found in the `probably approximately correct' (PAC)-Bayes approach to prediction in the machine learning literature (see alquier2021user, for a review), it is in the statistics and econometrics literature that this `loss-based prediction' has come to more formal maturity, including in terms of its theoretical validation. This includes Bayesian work on weighted combinations of predictions, such as, for example, Billio2013, casarin2015 , Pett2016, Bassetti2018, BASTURK2019, McAlinn2019{ and McAlinn2020,} where weights are updated via various predictive criteria, and the true model is not assumed to be one of the constituent models - i.e. an $\mathcal{M}$-open state of the world (Bernardo and Smith, 1994) is {implicitly }adopted. It also {includes a contribution} by loaiza2019focused, in which both single models and predictive mixtures are used to generate accurate Bayesian predictions in the presence of model misspecification; with both theoretical and numerical results highlighting the ability of the approach to out-perform conventional likelihood-based Bayesian predictive methods.
However, as a general rule, the existing approaches discussed above do not scale well to complex models with high-dimensional parameter spaces. In contrast, the current paper contributes to the above literature by providing a new approach for producing accurate loss-based predictions in high-dimensional problems. {We begin by defining a class of }flexible{\ predictive models{, }} conditional on a set of unknown parameters, {that are a plausible mechanism for generating probabilistic predictions. A prior distribution is placed over }the parameters of this predictive {class, and the prior then updated } to a posterior{\ via a criterion function that captures a { user-specified measure of predictive accuracy. That is, the conventional, and potentially misspecified likelihood-based update is eschewed, in favour of a function that is tailored to the predictive problem at hand; th}}e ultimate{{\ goal being to produce accurate predictions according to the measure that }}matters, without requiring knowledge of the true data generating mechanism.
In the spirit of the various generalized Bayesian inferential methods, in which likelihood functions are also replaced by alternative updating mechanisms (inter alia, chernozhukov2003mcmc, Zhang2006a{{, Zhang2006b, {jiang2008}, {bissiri2016general}}}{{{{, }giummole2017objective}}}, {{{GVI2019}}, }{miller2019robust, Syring2019} ,\ and pacchiardi2021generalized), we adopt a coherent update based on the exponential of a scaled sample loss; we refer to miller2021asymptotic for a {thorough discussion of} the large sample behavior of generalized Bayesian posteriors. As in loaiza2019focused the loss is, in turn, defined by a proper scoring rule ( gneiting2007probabilistic; gneiting2007strictly) that rewards a given form of predictive accuracy; for example, accurate prediction of extreme values. Given the high-dimensional nature of the resultant posterior, numerical treatment via `exact' Markov Chain Monte Carlo (MCMC) is computationally challenging; hence we adopt an `approximate' approach using variational principles. Since the posterior that results from an exponentiated loss was first denoted {as }a `Gibbs posterior' by Zhang2006a,\ we refer to the variational approximation of this posterior{\textbf{\ }}as the\textit{\ Gibbs variational posterior,} and the predictive distribution that results from this posterior, via the standard Bayesian calculus, as {the }\textit{Gibbs variational predictive } (hereafter, GVP). With a slight abuse of terminology, and when it is clear from the context, we also use the abbreviation GVP to reference the method of Gibbs variational prediction, or loss-based variational prediction \textit{per se.}
In an artificially simple `{toy' }example in which MCMC sampling of the exact Gibbs posterior {is feasible, we illustrate that there are negligible }differences between the out-of-sample results yielded by the GVP and those produced by the predictive based on MCMC sampling from the {Gibbs} posterior. We establish this result under both correct specification, in which the true data generating process matches the adopted predictive model, and under misspecification of the predictive model. The { correspondence} between the `approximate' and `exact' predictions in the correct specification case mimics that documented in frazier2019approximate, in which posterior approximations are produced by approximate Bayesian computation (ABC) and for the log score update (only).\ In the misspecified case, `strictly coherent' predictions are produced ( martin2020optimal), whereby a given GVP, constructed\ via the use of a particular scoring rule, is shown to perform best out-of-sample according to that same score when compared with a GVP constructed via some alternative scoring rule; with the numerical values of the average out-of-sample scores closely matching those produced by the exact Gibbs predictive. That is, building a generalized posterior via a given scoring rule yields superior predictive accuracy in that rule despite any inaccuracy in the measurement of posterior uncertainty that is induced by the variational approximation.
We then undertake more extensive Monte Carlo experiments to highlight the power of the approach in genuinely high-dimensional problems, with predictives based on: an autoregressive mixture model with 20 mixture components, and a neural network model used for illustration. An empirical analysis, in which GVP is used to produce accurate prediction intervals for the 4227 daily time series used in the M4 forecasting competition, illustrates the applicability of the method to reasonably large, and realistic data sets.
While the {numerical toy example} {demonstrates} that little predictive accuracy is lost when using the GVP, relative to the `exact' predictive, we {also }rigorously compare the theoretical behavior {of }the GVP and the potentially infeasible {exact} predictive. Specifically, we demonstrate that in large samples the GVP delivers predictions that are just as accurate as those obtained from the exact Gibbs predictive when measured according to the score used to construct the Gibbs posterior. We do this by proving that the GVP `merges' (in the sense of blackwell1962merging) with the exact Gibbs predictive. This merging result relies on novel posterior concentration results that extend existing concentration results for generalized posteriors (Alquier2016, alquier2020concentration, and yang2020alpha) to temporally dependent, potentially heavy-tailed data.
The remainder of the paper is structured as follows. In Section (ref) the loss-based paradigm for Bayesian prediction is defined, with its links to related segments of the literature (as flagged briefly above) detailed. In Section (ref) we detail how to construct the GVP, and provide theoretical verification of its accuracy. A numerical illustration of the ability of the variational method to yield essentially equivalent predictive accuracy to that produced via MCMC sampling is given in Section (ref) , using a low-dimensional example. The approach that we use in the implementation of GVP, including the choice of variational family, is briefly described, with further computational details of the stochastic gradient ascent (SGA) method used to perform the optimization are included in a supplementary appendix to the paper. We then proceed with the numerical illustration of the method in high-dimensional settings - using both artificially simulated and empirical data - in Sections (ref) and (ref) respectively; with the illustrations highlighting that, overall, the method `works' and `works well'. We conclude in Section (ref) with discussion of the implications of our findings, including for future research directions. Proof of all theoretical results are given in an appendix, and all computational details, and prior specifications are included in supplementary appendix to the paper.
Consider a stochastic process $\{Y_t:\Omega\rightarrow\mathcal{Y} , t\in \mathbb{N}\}$ defined on the complete probability space $(\Omega, \mathcal{F} ,P_0)$. Let $\mathcal{F}_t:=\sigma(Y_1,\dots,Y_t)$ denote the natural sigma-field, and let $P_0$ denote the infinite-dimensional distribution of the sequence $Y_1,Y_2,\dots$. Let $y_{1:n}=(y_1,\dots,y_n)^{\prime }$ denote a vector of realizations from the stochastic process.
Our goal is to use a particular collection of statistical models, adapted to $\mathcal{F}_{n}$, that describe the behavior of the observed data, to construct accurate predictions for the random variable $Y_{n+1}$. The parameters of the model are denoted by $\theta _{n}$, the parameter space by $\Theta _{n}\subseteq \mathbb{R}^{d_{n}}$, where the dimension $d_{n}$ could grow as $n\rightarrow \infty $ and $\Theta _{1}\subseteq \Theta _{2}\subseteq ...\subseteq \Theta _{n}$. For the notational simplicity, we drop\ the dependence of\ $\theta _{n}$ and $\Theta _{n}$ on $n$ in what follows. Let $\mathcal{P}^{(n)}$ be a generic class of one-step-ahead predictive models for $Y_{n+1}$, conditioned on the information $\mathcal{F} _{n}$ available at time $n$, such that $\mathcal{P}^{(n)}:=\{P_{\theta }^{(n)},\theta \in \Theta \}$.\footnote{{The treatment of scalar $Y_{t}$ and one-step-ahead prediction is for the purpose of illustration only, and all the methodology that follows can easily be extended to multivariate $Y_{t}$ and multi-step-ahead prediction in the usual manner.}} When $P_{\theta }^{(n)}(\cdot )$ admits a density with respect to the Lebesgue measure, we denote it by $p_{\theta }^{(n)}(\cdot )\equiv p_{\theta }(\cdot |\mathcal{F} _{n})$. The parameter ${\theta }$ thus\ indexes values in the predictive class, with ${\theta }$ taking values in the complete probability space $(\Theta ,\mathcal{T},\Pi )$, and where $\Pi $ {measures} our beliefs - {either prior or posterior - }about the unknown parameter ${\theta }$, and when {they} {exist} we denote {the respective densities} by $\pi (\theta )$ { and} $\pi (\theta |y_{1:n})$.
Denoting the likelihood function by $p_{\theta }(y_{1:n})$, the conventional approach to Bayesian prediction updates prior beliefs about ${\theta }$ via Bayes rule, to form the Bayesian posterior density,
in which we follow convention and abuse notation by writing this quantity as a density even though, strictly speaking, the density may not exist. The one-step-ahead predictive distribution is then constructed as
However, when the class of predictive models indexed by $\Theta $ does not contain the true predictive distribution there is no sense in which this conventional approach remains the `gold standard'. In such cases, the loss that underpins ((ref)) should be replaced by the particular predictive loss that matters for the problem at hand. That is, \ our prior beliefs about ${\theta }$ and, hence, about the elements $P_{\theta }^{(n)}$ in $\mathcal{P}^{(n)}$, need to be updated via a criterion function defined by a user-specified measure of predictive loss.
loaiza2019focused propose a method for producing Bayesian predictions using loss functions that specifically capture {the accuracy of} {density forecasts}. For $\mathcal{P}^{(n)}$ a convex class of {predictive distributions} on $(\Omega ,\mathcal{F})$, density prediction accuracy can be measured using the positively-oriented\ (i.e. higher is better) scoring rule $s:\mathcal{P}^{(n)}\times \mathcal{Y}\mapsto \mathbb{R} $, where the expected scoring rule under the true distribution $P_{0}$ is defined as
{Since} $\mathbb{S}(\cdot ,P_{0})$ is unattainable in practice, a sample estimate based on $y_{1:n}$ is used to define a sample criterion: for a given $\theta \in \Theta $, define sample average score as
Adopting the generalized updating rule proposed by bissiri2016general (see also giummole2017objective, holmes2017assigning, lyddon2019general, and Syring2019), loaiza2019focused distinguish elements in $ \Theta $ using
where the scale factor\ $w$ is obtained in a preliminary step using measures of predictive accuracy. This posterior explicitly weights elements of $\Theta $ according to their predictive accuracy in the scoring rule $s(\cdot ,\cdot )$. As such, the one-step-ahead generalized predictive,
will often outperform, in the chosen rule $s(\cdot ,\cdot )$, the likelihood (or {log-score})-based predictive $P_{\Pi }^{(n)}$ in cases where the model is misspecified. Given its explicit dependence {on the Gibbs posterior in ((ref))}, and as noted earlier, we refer to the predictive in ((ref)) as the (exact) Gibbs predictive.\footnote{We note that, in general the choice of $w$ can be $n$-dependent. However, we eschew this dependence to maintain notational brevity.}
While $\Pi _{w}(\cdot |y_{1:n})$ is our ideal posterior, it can be difficult to sample from if the dimension of $\theta $ is {even moderately} large, which occurs in situations with a large number of predictors, or in flexible models. Therefore, the exact predictive itself is not readily available in such cases. However, this does not invalidate the tenant on which $P_{\Pi _{w}}^{(n)}$ is constructed. Viewed in this way, we see that the problem of predictive inference via $P_{\Pi _{w}}^{(n)}$ could be solved if one were able to construct an accurate enough approximation to $\Pi _{w}(\cdot |y_{1:n})$. Herein, we propose to approximate $P_{\Pi _{w}}^{(n)}$ using variational Bayes (VB).
VB seeks to approximate $\Pi _{w}(\cdot |y_{1:n})$ by finding the closest member in a class of probability measures, denoted as $\mathcal{Q}$, to $\Pi _{w}(\cdot |y_{1:n})$ in a chosen divergence measure. The most common choice of divergence is the Kullback-Leibler (KL) divergence: for probability measures $ P,Q$, and $P$ absolutely continuous with respect to $Q$, the KL divergence is given by $\text{D} (P||Q)=\int \log (\mathrm{d} P/\mathrm{d} Q)\mathrm{d} P$. VB then attempts to produce a posterior approximation to $\Pi _{w}$ by choosing $Q$ to minimize:
In practice, VB is often conducted, in an equivalent manner, by maximizing the so-called evidence lower bound (ELBO). In the case where $\Pi_w$ and $Q$ admit densities $\pi_w$ and $q$, respectively, the ELBO is given by:
While other classes of divergences are used in variational inference, such as the R\'{e}nyi-divergence, the KL divergence is the most commonly encountered in the literature. Hence, in what follows we focus on this choice, but note here that other divergences can also be used.
The variational approximation to the Gibbs posterior, $P_{\Pi _{w}}^{(n)}$, can be defined as follows.
The GVP, $P_{Q}^{(n)}$, circumvents the need to construct $ P_{\Pi _{w}}^{(n)}$ via sampling the Gibbs posterior $\Pi _{w}(\cdot |y_{1:n}) $. In essence, we replace the sampling problem with an optimization problem for which reliable methods exist even if $\Theta $ is high-dimensional, and which in turn yields an approximation to $\Pi _{w}(\cdot |y_{1:n})$. Consequently, even though, as discussed earlier, it is not always feasible to {access} $P_{\Pi _{w}}^{(n)}$, via simulation from, $\Pi _{w}(\cdot |y_{1:n})$ in situations where $\Theta $ is high-dimensional, {access to the variational predictive }$P_{Q}^{(n)}$ remains feasible.
This variational approach to prediction is related to the `generalized variational inference'\ approach of GVI2019,\ but with two main differences. Firstly, our approach is focused on predictive accuracy, not on parameter inference. Our only goal is to produce density forecasts that are as accurate as possible in the chosen scoring rule $ s(\cdot ,y)$.\footnote{ Critically, since predictive accuracy is our goal, we are not concerned with the potential over/under-coverage of credible sets for the model unknowns built from generalized posteriors, or VB posteriors; see miller2021asymptotic for a theoretical discussion of posterior coverage in generalized Bayesian inference.} Secondly, our approach follows the idea of bissiri2016general and loaiza2019focused and targets the predictions built from the Gibbs posterior in (ref), rather than the exponentiated form of some general loss as in GVI2019. This latter point is certainly critical if inferential accuracy were still deemed to be important, since without the tempering constant $w$ that defines the posterior in (ref), the exponentiated\ loss function can be very flat and the posterior $\Pi _{w}(\cdot |y_{1:n})$ uninformative about $\theta $. \footnote{ See bissiri2016general, giummole2017objective, holmes2017assigning, lyddon2019general, Syring2019 and pacchiardi2021generalized{\ for various approaches to the setting of }$w$ { \ in inferential settings.}} It is our experience however ( loaiza2019focused), that in settings where predictive accuracy is the only goal, the choice of $w$ actually has little impact on the generation of accurate predictions via $P_{\Pi _{w}}^{(n)}$, as long as the sample size is reasonably large. Preliminary experimentation in the current setting, in which the variational predictive $P_{Q}^{(n)}$ is the target, suggests that this finding remains relevant, and as a consequence we have adopted a default value of $w=1$ in all numerical work. Further research is needed to deduce the precise impact of $w$ on $P_{Q}^{(n)}$, in particular for smaller sample sizes, and we leave this important topic for future work.\footnote{We note here that the work of wu2021calibrating attempts to calibrate the coverage of posterior predictives constructed from power posteriors. While the approach of wu2021calibrating is not directly applicable in our context, we conjecture that an appropriately modified version of their approach could be applied in our setting to derive posterior predictives with reasonable predictive coverage.}
By driving the updating mechanism by the {measure of predictive accuracy} that matters, it is hoped that {GVP produces} accurate predictions without requiring correct model specification. The following result demonstrates that, in large samples, the GVP, $P^{(n)}_Q$, produces predictions that are indistinguishable from the exact predictive, $P^{(n)}_{\Pi_{w}}$; thus, in terms of predictive accuracy, the use of a variational approximation to the Gibbs posterior has little impact on the accuracy of the posterior predictives. {For probability measures $P,Q$, let $d_{\text{TV} }\{P,Q\}$ denote the total variation distance. }
Theorem (ref) demonstrates that the difference between distributional predictions made using the GVP and the exact Gibbs predictive agree in the rule $s(\cdot ,\cdot ) $ in large samples. This type of result is colloquially known as a `merging' result (blackwell1962merging), and means that if {we take the exact Gibbs predictive as our benchmark for (Bayesian) prediction, then the GVP is as accurate as this possibly infeasible benchmark (at least in large samples).}
As a further result, we can demonstrate that GVP is as accurate as the infeasible frequentist `optimal [distributional] predictive' approach suggested in gneiting2007strictly. {Following gneiting2007strictly, define the `optimal' (frequentist) predictive within the class $\mathcal{P}^{(n)}$, and based on scoring rule $s(\cdot ,y)$, as $P_{\star }^{(n)}:=P(\cdot |\mathcal{F}_{n},\theta _{\star })$, where }
The following result demonstrates that, in large samples, predictions produced using {the GVP} are equivalent to those made using the optimal frequentist predictive $ P_{\star }^{(n)}$.
We remind the reader that, for the sake of brevity, statement of the assumptions, and proofs of all stated results, are given in Appendix (ref).
We now illustrate the behavior of GVP in a simple toy example for a financial return, $Y_{t}$, in which the exact Gibbs predictive is also accessible.\footnote{ We reiterate here, and without subsequent repetition, that all details of the prior specifications and the numerical steps required to produce the variational approximations, for this and the following numerical examples, are provided in {the supplementary appendices}.} This example serves {two purposes. First, it illustrates} {Theorem (ref)}, and, in so doing, highlights the benefits of GVP relative to prediction based on a misspecified likelihood function. {Secondly}, {the }GVP {is shown to }yield almost equal (average) out-of-sample predictions to those produced by the exact Gibbs posterior accessed via MCMC;{\ i.e. numerical support for Theorem (ref) is provided.}
The predictive class, $\mathcal{P}^{(n)}$, is defined by a generalized autoregressive conditional heteroscedastic GARCH(1,1) model with Gaussian errors, $Y_{t}=\theta _{1}+\sigma _{t}\varepsilon _{t},$\ $\varepsilon _{t} \overset{ i.i.d.}{\sim }N\left( 0,1\right) ,$\ $\sigma _{t}^{2}=\theta _{2}+\theta _{3}\left( Y_{t-1}-\theta _{1}\right) ^{2}+\theta _{4}\sigma _{t-1}^{2},$\ with $\theta =\left( \theta _{1},\theta _{2},\theta _{3},{ \theta }_{4}\right) ^{\prime }$. We adopt three alternative specifications for the true data generating process (DGP) : 1) a model that matches the Gaussian GARCH(1,1) predictive class; 2) a stochastic volatility model with leverage:
and 3) a stochastic volatility model with a smooth transition in the volatility autoregression:
where $g(x)=\left( 1+\exp \left( -2x\right) \right) ^{-1}$. DGP 1) defines a \ correct specification\ setting, whilst DGPs 2) and 3) characterize different forms of\ misspecification.
Denoting the predictive density function associated with the Gaussian GARCH (1,1) model, evaluated at the observed\ $y_{t+1}$,\ as $ p_{\theta }(y_{t+1}|\mathcal{F}_{t})$, we implement GVP using four alternative forms of (postively-oriented) scoring rules:
where $l_{t+1}\ $and $u_{t+1}$\ denote the $100(\frac{\alpha }{2})\%$\ and $ 100(1-\frac{\alpha }{2})\%\ $predictive quantiles. The log-score (LS) in ( (ref)) is a `local' scoring rule, attaining a high value if the observed value, $y_{t+1}$, is in the high density region of $p_{\theta }(y_{t+1}|F_{t})$. The continuously ranked probability score (CRPS) in ((ref)) (gneiting2007strictly) is, in contrast, sensitive to distance,\ and rewards the assignment of high predictive mass near to the realized $y_{t+1}$, rather than just at that value. The score in ((ref)) is the censored likelihood score (CLS)\ of diks2011likelihood, which rewards predictive accuracy over any pre-specified region of interest $A$\ ($ A^{c}$\ denoting the complement). We use the score for $A$\ defining the lower and upper tails of the predictive distribution, as determined in turn by the 10%, 20%, 80% and 90% percentiles\ of the empirical distribution of $Y_{t}$, labelling these cases hereafter as CLS$10$, CLS$20,$ CLS$80$ and CLS$90.$ A high value of this score in any particular instance thus reflects a predictive distribution that accurately predicts extreme values of the financial return. The last score considered is the{\ interval score (IS)} in ((ref)), which is designed to measure the accuracy of the $100(1-\alpha )\%$\ predictive interval where $ \alpha =0.05$. This score rewards narrow intervals with accurate coverage. All components of ((ref)) have closed-form solutions for the (conditionally) Gaussian predictive model, as does the integral in ((ref)) and the bounds in ((ref)).
In total then, seven distinct scoring rules are used to define the sample criterion function in ((ref)). In what follows we reference these seven criterion functions, and the associated Gibbs posteriors using the notation\ $S_{n}^{j}\left( \theta \right) =\sum_{t=0}^{n-1}s^{j}\left( P_{\theta }^{(t)},y_{t+1}\right) $,\ for\ $j=\{\text{LS, CRPS, CLS10, CLS20, CLS80, }${CLS90, IS}$ \}$.
Given the {low dimension} of the predictive model it is straightforward to use an MCMC scheme to sample from the exact Gibbs posterior $\pi _{w}^{j}({ \theta }|y_{1:n})\propto \exp \{wS_{n}^{j}(\theta )\}\pi \left( \theta \right) , $ where $\pi _{w}^{j}(\theta |y_{1:n})$ is the exact Gibbs posterior density in ((ref)) computed under scoring rule $j.$ As noted earlier, in this and all following {numerical} work we set $w=1$. For each of the $j$ posteriors, we initialize the chains using a burn-in period of $ 20000$ periods, and retain the next $M=20000$\ draws $\theta ^{(m,j)}\sim \pi _{w}^{j}(\theta |y_{1:n}),$ $m=1,\dots ,M$. The posterior draws are then used to estimate the exact Gibbs predictive in ((ref)) via $\hat{P} _{\Pi _{w}}^{(n,j)}=\frac{1}{M}\sum_{m=1}^{M}P_{\theta ^{(m,j)}}^{(n)}. $
To perform GVP, we first need to produce the variational approximation of $ \pi _{w}^{j}(\theta |y_{1:n}).$ This is achieved in several steps. First, the parameters of the GARCH(1,1) model are transformed to the real line. With some abuse of notation, we re-define here the GARCH(1,1) parameters introduced in the previous section with the superscript $r$\ to signal `raw'. The parameter vector $\theta $ is then re-defined as $\theta =\left( \theta _{1},\theta _{2},\theta _{3},\theta _{4}\right) ^{\prime }=\left( \theta _{1}^{r},\log (\theta _{2}^{r}),\Phi _{1}^{-1}\left( \theta _{3}^{r}\right) ,\Phi _{1}^{-1}\left( \theta _{4}^{r}\right) \right) ^{\prime }$, where $\Phi _{1}$\ denotes the {(univariate)} normal cumulative distribution function (cdf). The next step involves approximating $\pi _{w}^{j}(\theta |y_{1:n})$ for the re-defined $\theta $. We adopt the mean-field variational family (see for example blei2017variational ), with a product-form Gaussian density $q_{\lambda }\left( \theta \right) =\prod_{i=1}^{4}\phi _{1}\left( {\theta }_{i};\mu _{i},d_{i}^{2}\right) $, where $\lambda =\left( \mu ^{\prime },d^{\prime }\right) ^{\prime }$ is the vector of variational parameters, {comprised of }mean and variance vectors $ \mu =\left( \mu _{1},\mu _{2},\mu _{3},\mu _{4}\right) ^{\prime }$ and $ d=\left( d_{1},\dots ,d_{4}\right) ^{\prime }$\ respectively, {and }$\phi _{1}$ {denotes the {(univariate) }normal probability density function (pdf). }Denote as $Q_{\lambda }$ and $\Pi _{w}^{j}$ the {distribution functions associated with} $q_{\lambda }\left( \theta \right) $ and $\pi _{w}^{j}(\theta |y_{1:n})$, respectively. The approximation is then calibrated by solving the maximization problem
where, $\mathcal{L}(\lambda ):=\text{ELBO}\left[ Q_{\lambda }||\Pi _{w}^{j} \right] $. Remembering the use of the notation $S_{n}^{j}\left( \theta \right) $ to denote ((ref)) for scoring rule $j$, ((ref)) (for case $j$) becomes $\mathcal{L}(\lambda )=\mathbb{E}_{q_{\lambda }}\left[ wS_{n}^{j}\left( {\theta }\right) +\log \pi \left( {\theta }\right) -\log q_{\lambda }\left( {\theta }\right) \right] .$ Optimization is performed via SGA, as described in Section (ref) of the supplementary appendix. Once calibration of $\tilde{\lambda}$ is completed, the GVP is estimated as $\hat{P}_{Q}^{(n,j)}=\frac{1}{M}\sum_{m=1}^{M}P_{\theta ^{(m,j)}}^{(n)}$ with $\theta ^{(m,j)}\overset{i.i.d.}{\sim }q_{\tilde{ \lambda}}\left( {\theta }\right) $. To calibrate $\tilde{\lambda}$ $10000$ VB iterations are used, and to estimate the variational predictive we set $ M=1000$.
To produce the numerical results we generate a times series of length $T=6000 $ from the true DGP. Then, we {perform} an expanding window exercise from $n=1000$ to $n=5999$. For $j\in \{\text{LS, CRPS, CLS10, CLS20, CLS80, CLS90, IS\}}$ we construct the predictive densities $\hat{P}_{\Pi _{w}}^{(n,j)} $ and $\hat{P}_{Q}^{(n,j)}$ as outlined\ above. Then, for $i\in \{\text{LS, CRPS, CLS10, CLS20, CLS80, CLS90, IS\}}$ we compute the measures of out-of-sample predictive accuracy $S_{\Pi _{w}}^{i,j,n}=s^{i}\left( \hat{P} _{\Pi _{w}}^{(n,j)},y_{n+1}\right) $ and $S_{Q}^{i,j,n}=s^{i}\left( \hat{P} _{Q}^{(n,j)},y_{n+1}\right) $. Finally, we compute the average out-of-sample scores $S_{\Pi _{w}}^{i,j}=\frac{1}{5000}\sum_{n=1000}^{5999}S_{\Pi _{w}}^{i,j,n}$ and $S_{Q}^{i,j}=\frac{1}{5000} \sum_{n=1000}^{5999}S_{Q}^{i,j,n}$. The results are tabulated and discussed in the following section.
The results of the simulation exercise are recorded in Table (ref). Panels A to C record the results for Scenarios 1) to 3), with the average out-of-sample scores associated with the exact Gibbs predictive (estimated via MCMC) appearing in the left-hand-side panel and the average scores for GVP appearing in the right-hand-side panel. All values on the diagonal of each sub-panel correspond to the case where $i=j$. In the misspecified case, numerical validation of the asymptotic result that the GVP concentrates onto the optimal predictive ({Theorem (ref)}) occurs if the largest values (bolded) in a column appear on the diagonal (`strict coherence' in the language of martin2020optimal). In the correctly specified case, in which all proper scoring rules will, for a large enough sample, pick up {\ the one} true model (gneiting2007strictly), we would expect all values in given column to be very similar to one another, differences reflecting sampling variation only. Finally, validation of the theoretical property of merging between the exact and variational Gibbs predictive ({Theorem (ref)}) occurs if the corresponding results in all left- and right-hand-side panels are equivalent.
As is clear, there is almost uniform numerical validation of {both theoretical results}. Under misspecification Scenario 2), (Panel B) the GVP results are strictly coherent (i.e. all bold values lie on the diagonal), with the exact Gibbs predictive results equivalent to the corresponding GVP results to three or four decimal places. The same broad findings obtain under misspecification Scenario 3), apart from the fact that the {IS updates} are second best (to the log-score; and then only just) in terms of the out-of-sample {IS measure}. In Panel A on the other hand, we see the expected (virtual) equivalence of all results in a given column, reflecting the fact that all Gibbs predictives (however estimated) are concentrating onto the true predictive model and, hence, have identical out-of-sample performance. Of course, for a large but still finite number of observations, we would expect the log-score to perform best, due to the efficiency of the implicit maximum likelihood estimator underpinning the results and, to all intents and purposes this is exactly what we observe in Panels A.1 and A.2.
In summary, GVP performs as anticipated, and reaps distinct benefits in terms of predictive accuracy. Any inaccuracy in the measurement of parameter uncertainty also has negligible impact on the finite sample performance of GVP relative to an exact comparator. {In Section (ref)}, we extend the investigation into design settings that mimic the high-dimensional problems to which we would apply the variational approach in practice, followed by an empirical application - again using a high-dimensional predictive model - in Section (ref).
In this section we demonstrate the application of GVP in two realistic simulated examples. In both cases the assumed predictive model is high-dimensional and the exact Gibbs posterior, even if accessible in principle via MCMC, is challenging from a computational point of view. The mean-field class is adopted in both cases to produce variational approximations to the Gibbs posterior. The simulation design for each example (including the choice of $w=1$) mimics that for the toy example, apart from the obvious changes {made to} the true DGP and the assumed predictive model, some changes {in} the size of the estimation and evaluation periods, plus the absence of comparative exact results. For reasons of computational burden we remove the CRPS update from consideration in the first example.
In this example we adopt a true DGP in which $Y_{t}$ evolves according to the logistic smooth transition autoregressive (LSTAR) process proposed in terasvirta1994specification:
where $\varepsilon _{t}\overset{i.i.d.}{\sim}t_{\nu }$, and $t_{\nu }$ denotes the standardised Student-t distribution with $\nu $ degrees of freedom. This model has the ability to produce data that exhibits a range of complex features. For example, it not only allows for skewness in the marginal density of $Y_{t}$, but can also produce temporal dependence structures that are asymmetric. We thus use this model as an illustration of a complex DGP whose characteristics are hard to replicate with simple parsimoneous models. Hence the need to adopt a highly parameterized predictive model; plus the need to acknowledge that even that representation will be misspecified.
The assumed predictive model is based on the flexible Bayesian non-parametric structure proposed in antoniano2016nonparametric. The predictive distribution for $Y_{t}$, conditional on the observed $y_{t-1}$, is constructed from a mixture of $K=20$\ Gaussian autoregressive (AR) models of order one as follows:
with time-varying mixture\ weights
where $\mu _{k}=\frac{\beta _{k,0}}{1-\beta _{k,1}}$ and $s_{k}^{2}=\frac{ \sigma _{k}^{2}}{1-\beta _{k,1}^{2}}$. Denoting $\tau=\left( \tau_{1},\dots ,\tau_{K}\right) ^{\prime }$, $\beta _{0}=\left( \beta _{1,0},\dots ,\beta _{K,0}\right) ^{\prime }$, $\beta _{1}=\left( \beta _{1,1},\dots ,\beta _{K,1}\right) ^{\prime }$ and $\sigma =\left( \sigma _{1},\dots ,\sigma _{K}\right) ^{\prime }$, then the full vector of unknown\ parameters is $\theta =\left( \mu ,\tau^{\prime },\beta _{0}^{^{\prime}},\beta _{1}^{\prime },\right.$ $\left.\sigma ^{\prime }\right) ^{\prime }$, which comprises $1+\left( 4\times 20\right) =81$ elements. GVP is a natural and convenient alternative to exact Gibbs prediction in this case.
In the second example we consider a true DGP in which the dependent variable $Y_{t}$ has a complex non-linear relationship with a set of covariates. Specifically, the time series process $\{Y_{t}\}_{t=1}^{T}$, is determined by a three-dimensional stochastic process $\{X_{t}\}_{t=1}^{T}$, with $ X_{t}=\left( X_{1,t},X_{2,t},X_{3,t}\right) ^{\prime }$. The first two covariates are jointly distributed as $(X_{1,t},X_{2,t})^{\prime }\overset{ i.i.d.}{\sim }N\left( 0,\Sigma \right) $. The third covariate $X_{3,t}$, independent of the former two, is distributed according to an AR(4) process so that $X_{3,t}=\sum_{i=1}^{4}\alpha _{i}X_{3,t-i}+\sigma \varepsilon _{t}$ , with $\varepsilon _{t}\overset{i.i.d.}{\sim }N\left( 0,1\right) $. The variable $Y_{t}$ is then given by $Y_{t}=X_{t}^{\prime }\beta _{t}$, where $ \beta _{t}=\left( \beta _{1,t},\beta _{2,t},\beta _{3,t}\right) ^{\prime }$ is a three dimensional vector of time-varying coefficients, with $\beta _{i,t}=b_{i}+a_{i}F\left( X_{3,t}\right) $, and $F$ denotes the marginal distribution function of $X_{3,t}$ induced by the AR(4) model.
This particular choice of DGP has two advantages. First, given the complex dependence structure of the DGP (i.e. non-linear cross-sectional dependence as well as temporal dependence), it would be difficult, once again, to find a simple parsimoneous predictive model that would adequately capture this structure; hence motivating the need to adopt a high-dimensional model and to resort to GVP. Second, because it includes several covariates, we can assess how GVP performs for varying informations sets. For example, we can establish if the overall performance of GVP improves with an expansion of the conditioning information set, as we would anticipate; and if the inclusion of a more complete conditioning set affects the occurrence of strict coherence, or not.
A flexible model that has the potential to at least partially capture several features of the true DGP\ is a Bayesian feed-forward neural network that takes $Y_{t}$ as the dependent variable and some of its lags, along with observed values of\ $X_{t}$, $x_{t}=\left( x_{1,t},x_{2,t},x_{3,t}\right) ^{\prime },$ as the vector of\ independent inputs, $z_{t}$ (see for instance hernandez2015probabilistic). The predictive distribution is defined as $P_{\theta }^{(t-1)}=\Phi _{1}\left( Y_{t};g\left( z_{t};\omega \right) ,\sigma _{y}^{2}\right) .$ The mean function $g\left( z_{t};\omega \right) $ denotes a feed-forward neural network with $q=2$ layers, $r=3$ nodes in each layer, a $p$-dimensional input vector $z_{t}$, and a $d$ -dimensional parameter vector $\omega $ with $d=r^{2}(q-1)+r(p+q+1)+1$. It allows for a flexible non-linear relationship between $Y_{t}$ and the vector of observed covariates $z_{t}$. Defining $c=\log (\sigma _{y})$, the full parameter vector of this predictive class, $\theta =\left( \omega ^{\prime },c\right) ^{\prime }$, is of dimension $d+1$.
A time series of length $T=2500$ for $Y_{t}$ is generated from the LSTAR model in (ref), with the true parameters set as: $\rho _{1}=0$, $\rho _{2}=0.9$, $\gamma =5$, $c=0$, $\sigma _{\varepsilon }=1$ and $\nu =3$ . With exception of the CRPS, which cannot be evaluated in closed-form for the mixture predictive class and is not included in the exercise as a consequence, the same scores in the toy example are considered. With reference to the simulation steps given earlier, the initial estimation window is $500$ observations; hence the predictive accuracy results are based on an evaluation sample of size $2000.$
The results in Table (ref) are clear-cut. With one exception (that of CLS10), the best out-of-sample performance, according to a given measure, is produced by the version of GVP based on that same scoring rule. That is, the GVP results are almost uniformly strictly coherent: a matching of the update rule with the out-of-sample measure produces the best predictive accuracy in that measure, almost always.
In this case, we generate a time series of length $T=4000$ from the model discussed in Section (ref), with specifications: $\Sigma _{11}=1$, $ \Sigma _{22}=1.25$, $\Sigma _{12}=\Sigma _{21}=0.5$, $\sigma ^{2}=0.2$, $ \alpha _{1}=0.5$, $\alpha _{2}=0.2$, $\alpha _{3}=0.15$, $\alpha _{4}=0.1$, $ a_{1}=1.3$, $b_{1}=0$, $a_{2}=-2.6$, $b_{2}=1.3$, $a_{3}=-1.5$ and $ b_{3}=1.5 $. These settings generate data with a negatively skewed empirical distribution, a non-linear relationship between the observations on $Y_{t}$ and $X_{3,t}$, and autoregressive behavior in $Y_{t}$.
To assess if the performance of GVP is affected by varying information sets, we consider four alternative specifications for the input vector $z_{t}$ in the assumed predictive model. These four specifications (labelled as Model 1, Model 2, Model 3 and Model 4, respectively) are: $z_{t}=y_{t-1}$, $ z_{t}=\left( y_{t-1},x_{1,t}\right) ^{\prime }$, $z_{t}=\left( y_{t-1},x_{2,t}\right) ^{\prime }$ and $z_{t}=\left( y_{t-1},x_{1,t},x_{2,t}\right) ^{\prime }$. The dimension of the parameter vector for each of these model specifications is $d+1=23$, $d+1=26$, $d+1=26$ and $d+1=29$, respectively. Given that the assumed predictive class is Gaussian, all scoring rules used in the seven updates, and for all four model specifications, can be evaluated in closed form. Referencing the simulation steps given earlier, the initial estimation window is $2000$ observations; hence the predictive accuracy results are based on an evaluation sample of size $2000$ as in the previous example.
With reference to the results in Table (ref) we can make two observations. First, as anticipated, an increase in the information set produces higher average scores out-of-sample, for all updates; with the corresponding values increasing steadily and uniformly as one moves from the results in Panel A (based on $z_{t}=y_{t-1}$) to those in Panel D (based on the largest information, $z_{t}=\left( y_{t-1},x_{1,t},x_{2,t}\right) ^{\prime }$) set. Secondly however, despite the improved performance of all versions of GVP as the information set is increased to better match that used in the true DGP, strict coherence still prevails. That is, `focusing' on the measure that matters in the construction of the GVP update, still reaps benefits, despite the reduction in misspecification.
The Makridakis 4 (M4) forecasting competition was a forecasting event first organised by the University of Nicosia and the New York University Tandon School of Engineering in 2018. The competition sought submissions of point and interval predictions at different time horizons, for a total of 100,000 times series of varying frequencies. The winner of the competition in a particular category (i.e. point\ prediction or interval prediction) was the submission that{\ achieved the best average out-of-sample predictive accuracy according to the measure of accuracy that defined that category, over all horizons and all series.\footnote{ Details of all aspects of the competition can be found via the following link:\newline https://www.m4.unic.ac.cy/wp-content/uploads/2018/03/M4-Competitors-Guide.pdf }}
We gauge the success of our method of\ distributional prediction in terms of the measure used to rank the interval forecasts in the competition, namely the IS. We focus on one-step-ahead prediction of each of the 4227 time series of daily frequency. Each of these series is denoted by $\left\{ Y_{i,t}\right\} $ where $i=1,\dots ,4227$ and $ t=1,\dots ,n_{i}$, and the task is to construct a predictive interval for $ Y_{i,n_{i}+1}$ based on GVP. We adopt the mixture distribution in (ref) as the predictive model, being suitable as it is to capture the stylized features of high frequency data. However, the model is only appropriate for stationary time series and most of the daily time series exhibit non-stationary patterns. To account for this, we model the differenced series $\left\{ Z_{i,t}=\Delta ^{d_{i}}Y_{i,t}\right\} $, where $ d_{i}$ indicates the differencing order. The integer $d_{i}$ is selected by sequentially applying the KPSS test (KPSS) to $\Delta ^{j}y_{i,t}$ for $j=0,\dots ,d_{i}$, with $d_{i}$ being the first difference at which the null hypothesis of no unit root is not rejected.\footnote{ To do this we employ the autoarima function in the forecast R package.} To construct the predictive distribution of $Z_{i,n_{i}+1}$ we first obtain $ M=5000$ draws $\{\theta _{i}^{(m)}\}_{m=1}^{M}$ from the Gibbs variational posterior $q_{\tilde{\lambda}}\left( \theta _{i}\right) $ that is based on the IS.\footnote{ Once again, a value of $w=1$ is used in defining the Gibbs posterior and, hence, in producing the posterior approximation.} Draws $ \{z_{i,n_{i}+1}^{(m)}\}_{m=1}^{M}$ are then obtained from the corresponding predictive distributions $\{P_{\theta _{i}^{(m)}}^{(n_{i}+1)}\}_{m=1}^{M}$. The draws of $Z_{i,n_{i}+1}$ are then transformed into draws of $ Y_{i,n_{i}+1}$ and a predictive distribution for $Y_{i,n_{i}+1}$ produced using kernel density estimation; assessment is performed in terms of the accuracy of the prediction interval for $Y_{i,n_{i}+1}.$ To account for different units of measurement, the IS of each series is then divided by the constant $\frac{1}{(n_i-1)}\sum_{t=2}^{n_i}|Y_{i,t}-Y_{i,t-1}|$.
Table (ref) documents the predictive performance of our approach relative to the competing methods makridakis2020m4 for details on these methods). The first column corresponds to the average IS. In terms of this measure the winner is the method proposed by smyl2019hybrid, while our method is ranked 11 out of 12. The ranking is\ six out of 12 in terms of the median IS, recorded in column two. However, it is worth remembering that our method aims to achieve high individual, and not aggregate (or average) predictive interval accuracy across series. A more appropriate way of judging the effectiveness of our approach is thus to count the number of series for which each method performs best. This number is reported in column three of the table. In terms of this measure, with a total number of 858 series, our approach is only outperformed by the method proposed in doornik2019card. In contrast, although smyl2019hybrid is the best method in terms of mean IS, it is only best at predicting 130 of the series in the dataset. It is also important to mention that most of the competing approaches have more complex assumed predictive classes than our mixture model. Even the simpler predictive classes like the ARIMA and ETS models have model selection steps that allow them to capture richer Markovian processes than the mixture model which, by construction, is based on an autoregressive process of order one.
In summary, despite a single specification being defined for the predictive class for all 4227 daily series, the process of driving the Bayesian update via the IS rule has still yielded the best\ predictive results in a very high number of cases, with GVP beaten in this sense by only one other competitor. Whilst the predictive model clearly matters, designing a bespoke updating mechanism to produce the predictive distribution is shown to still reap substantial benefits.
We have developed a new approach for conducting loss-based Bayesian prediction in high-dimensional models. Based on a variational approximation to the Gibbs posterior defined, in turn, by the predictive loss that is germane to the problem at hand, the method is shown to produce predictions that minimize that loss out-of-sample. `Loss' is characterized in this paper by {positively-oriented proper scoring} rules designed to reward the accuracy of predictive probability density functions for a continuous random variable. Hence, loss minimization translates to maximization of an expected score. However, in principle, any loss function in which predictive accuracy plays a role could be used to define the Gibbs posterior. Critically, we have proven theoretically, and illustrated numerically, that for a large enough sample there is no loss incurred in predictive accuracy as a result of approximating the Gibbs posterior.
{In comparison with the standard approach based on a likelihood-based Bayesian posterior, }our Gibbs variational predictive approach is ultimately aimed at generating accurate predictions in the realistic empirical setting where the predictive model and, hence, the likelihood function is misspecified. Gibbs variational prediction enables the investigator to break free from the shackles of likelihood-based prediction, {and to drive } predictive outcomes according to the form of predictive accuracy that matters for the problem at hand; and all with theoretical validity guaranteed.
We have focused in the paper on particular examples where the model used to construct the Gibbs variational predictive is observation-driven. Extensions to parameter-driven models (i.e., hidden Markov models, or state space models) require different approaches with regard to both the implementation of variational Bayes, and in establishing the asymptotic properties of the resulting posteriors and predictives. This is currently the subject of on-going work by the authors, frazier2021note, and we reference tran2017variational and quiroz2018gaussian for additional discussion of this problem.
This paper develops the theory of Gibbs variational prediction for classes of models that are general and flexible enough to accommodate a wide range of data generating processes, and which, under the chosen loss function, are smooth enough to permit a quadratic expansion. This latter condition restricts the classes of models, and loss functions, under which our results are applicable. For example, the discrete class of models studied in douc2013ergodicity may not be smooth enough to deduce the validity of such an expansion.
While our approach to prediction focuses on accuracy in a given forecasting rule, which is related to the notion of sharpness in the probabilistic forecasting literature (see, gneiting2007probabilistic), it does not pay attention to the calibration of such predictive densities. That is, there is no guarantee that the resulting predictive densities have valid frequentist coverage properties. If calibration is a desired property, it is possible to augment our approach to prediction with the recently proposed approach by wu2021calibrating, which calibrates generalized predictives. The amalgamation of these two procedures should produce accurate predictions that are well-calibrated in large samples. We leave an exploration of this possibility to future research.
\acks{We would like to thank various participants at: the `ABC in Svalbard' Workshop, April 2021, the International Symposium of Forecasting, the Alan Turing Institute, and RIKEN reading group, June 2021, the European Symposium of Bayesian Econometrics, September 2021, and the Vienna University of Business, December 2021, for helpful comments on earlier drafts of the paper. This research has been supported by Australian Research Council (ARC) Discovery Grants DP170100729 and DP200101414. Frazier was also supported by ARC Early Career Researcher Award DE200101070.}