Loss-Based Variational Bayes Prediction

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\}$.¹¹1The 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,

$$\pi(\theta|y_{1:n})=\frac{p_{\theta}\left(y_{1:n}\right)\pi(\theta)}{\int_{\Theta}p_{\theta}\left(y_{1:n}\right)\pi(\theta)d\theta},$$

(1)

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

$$P_{\Pi}^{(n)}:=\int_{\Theta}P_{\theta}^{(n)}\pi(\theta|y_{1:n})\mathrm{d}\theta.$$

(2)

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 (2) 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.

Loaiza-Maya et al. (2020) 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

$$\mathbb{S}(\cdot,P_{0}):=\int_{y\in\Omega}s(\cdot,y)\mathrm{d}P_{0}(y).$$

(3)

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

$$S_{n}(\theta):=\sum_{t=0}^{n-1}s(P_{\theta}^{(t)},y_{t+1}).$$

(4)

Adopting the generalized updating rule proposed by Bissiri et al. (2016) (see also Giummolè et al., 2017, Holmes and Walker, 2017, Lyddon et al., 2019, and Syring and Martin, 2019), Loaiza-Maya et al. (2020) distinguish elements in $\Theta$ using

$$\pi_{w}(\theta|y_{1:n})=\frac{\exp\left[wS_{n}\left(\theta\right)\right]\pi(\theta)}{\int_{\Theta}\exp\left[wS_{n}\left(\theta\right)\right]\pi(\theta)\mathrm{d}\theta},$$

(5)

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,

$$P_{\Pi_{w}}^{(n)}:=\int_{\Theta}P_{\theta}^{(n)}\pi_{w}(\theta|y_{1:n})\mathrm{d}\theta,$$

(6)

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 (5), and as noted earlier, we refer to the predictive in (6) as the (exact) Gibbs predictive.²²2We 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:

$$\text{D}\left(Q||\Pi_{w}\right)=\int\log\left[{\mathrm{d}Q}/{\mathrm{d}\Pi_{w}(\cdot|y_{1:n})}\right]\mathrm{d}Q.$$

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:

$$\text{ELBO}[Q||\Pi_{w}]:=\mathbb{E}_{q}[\log\left\{\exp\left[wS_{n}(\theta)\right]\pi(\theta)\right\}]-\mathbb{E}_{q}[\log\{q(\theta)\}].$$

(7)

While other classes of divergences are used in variational inference, such as the Ré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 Knoblauch et al. (2019), 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)$.³³3Critically, 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 Miller (2021) for a theoretical discussion of posterior coverage in generalized Bayesian inference. Secondly, our approach follows the idea of Bissiri et al. (2016) and Loaiza-Maya et al. (2020) and targets the predictions built from the Gibbs posterior in (5), rather than the exponentiated form of some general loss as in Knoblauch et al. (2019). 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 (5), the exponentiated loss function can be very flat and the posterior $\Pi_{w}(\cdot|y_{1:n})$ uninformative about $\theta$.⁴⁴4See Bissiri et al. (2016), Giummolè et al. (2017), Holmes and Walker (2017), Lyddon et al. (2019), Syring and Martin (2019) and Pacchiardi and Dutta (2021) for various approaches to the setting of $w$  in inferential settings. It is our experience however (Loaiza-Maya et al., 2020), 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.⁵⁵5We note here that the work of Wu and Martin (2021) attempts to calibrate the coverage of posterior predictives constructed from power posteriors. While the approach of Wu and Martin (2021) 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 1 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 (Blackwell and Dubins, 1962), 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 Gneiting and Raftery (2007). Following Gneiting and Raftery (2007), 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

$$\theta_{\star}:=\arg\max_{\theta\in\Theta}\mathcal{S}(\theta),\text{ {with} }\mathcal{S}(\theta)=\operatornamewithlimits{lim\,}_{n\rightarrow\infty}\mathbb{E}\left[S_{n}(\theta)/n\right].$$

(8)

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 A.

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 (5) 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 (6) 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 Blei et al., 2017), 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

$$\tilde{{\lambda}}=\text{ }\underset{\lambda\in\Lambda}{\arg\max}\ \mathcal{L}(\lambda),$$

(13)

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 (4) for scoring rule $j$, (7) (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 S2 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 1. 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 2) occurs if the largest values (bolded) in a column appear on the diagonal (‘strict coherence’ in the language of Martin et al., 2022). In the correctly specified case, in which all proper scoring rules will, for a large enough sample, pick up  the one true model (Gneiting and Raftery, 2007), 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 1) 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 5, 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 6.

In this example we adopt a true DGP in which $Y_{t}$ evolves according to the logistic smooth transition autoregressive (LSTAR) process proposed in Teräsvirta (1994):

$$Y_{t}=\rho_{1}Y_{t-1}+\rho_{2}\left\{\frac{1}{1+\exp\left[-\gamma\left(Y_{t-1}-c\right)\right]}\right\}y_{t-1}+\sigma_{\varepsilon}\varepsilon_{t},$$

(14)

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 Antoniano-Villalobos and Walker (2016). 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:

$$P_{\theta}^{(t-1)}=\sum_{k=1}^{K}\tau_{k,t}\Phi_{1}\left[Y_{t}-\mu;\beta_{k,0}+\beta_{k,1}(y_{t-1}-\mu),\sigma_{k}^{2}\right],$$

(15)

with time-varying mixture weights

$$\tau_{k,t}=\frac{\tau_{k}\phi_{1}\left(y_{t-1}-\mu;\mu_{k},s_{k}^{2}\right)}{\sum_{j=1}^{K}\tau_{j}\phi_{1}\left(y_{t-1}-\mu;\mu_{j},s_{j}^{2}\right)},$$

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 Hernández-Lobato and Adams, 2015). 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$.

For some positive sequence $\epsilon_{n}\rightarrow 0$, define $\Theta(\epsilon_{n}):=\{\theta\in\Theta:d(\theta,\theta_{\star})\leq\epsilon_{n}\}$ as an $n$-dependent neighbourhood of $\theta_{\star}$ and let $\Theta^{c}(\epsilon_{n})$ denote its complement.