5.3 COMPUTATIONAL METHODS

5.3 COMPUTATIONAL METHODS

Obtaining a Bayesian point forecast is straightforward in cases, such as those examined above, where we have both conjugate priors and point forecasts that depend on a known function of the parameter distribution. In such cases it is also simple to calculate predictive distributions or density forecasts. Such distributions can be useful in many ways—we can compute averages of posterior distributions of the model parameters, compute the average loss to understand the properties of the forecast, and construct density forecasts or make probability statements on how likely it is that the outcome falls in a specified range.

More generally, however, numerical methods are required to make these calculations. Numerical methods can be employed to obtain estimates of posterior distributions for the parameters $\pi (\theta |z)$ , predictive densities $P^P(y|z)$ , as well as summary statistics of these distributions. For example, the posterior density of the model parameters can be used to compute the mean outcome and construct onestep-ahead forecasts under MSE loss. Alternatively, $\pi (\theta |z)$ could be used to generate the predictive density, the mean of which can be used in multistep forecasting.

The standard problem in Bayesian estimation is to compute the posterior density of the parameters, $\pi (\theta |z)$ . This is frequently very involved with different methods tailored to the particular model that needs to be evaluated. Since many texts are devoted to this subject, we give only a cursory overview of the methods to help the reader access the literature. Later chapters are clearer on the specific methods used for a range of forecasting problems.

The basic idea is to design a computational method for constructing draws from the posterior distribution π $\left(\theta |z\right)$ even though its exact density is unknown in cases where we are unable to describe it as a known distribution of the parameters of the prior and the data. Instead of directly working with the posterior, typically $p(z|\theta )\pi (\theta )$ is examined without normalizing by the marginal distribution of z. Normalization of this density is easily achieved using numerical approximation. In any particular application, $p(z|\theta )$ and $\pi (\theta )$ are both known since the forecaster must decide on a likelihood for the observed data z and a prior over the parameters, θ.

The next step is to sample from the distribution $p(z|\theta )\pi (\theta )$ , given the data. This involves sampling from the prior multiplied by the likelihood, both of which are fully known as z is observed and θ is a draw from the prior. This can then be normalized using the draws. To enhance efficiency, often steps are taken to narrow down the support of the posterior so that sampling is undertaken in the correct range. Unfortunately, this approach is not very efficient or easy to use when $\theta$ is large-dimensional and contains more than a few parameters. For this common case, numerical methods that draw from a distribution that converges to the posterior are employed, rather than drawing θ directly from the posterior. Monte Carlo

Markov Chain (MCMC) methods typically involve refinements such as the Gibbs sampler or the Metropolis algorithm for drawing from the underlying distributions.

The Gibbs sampler is useful for situations where we cannot find a closed form for $\pi (\theta |z)$ but can split θ into subsets, say $\theta =({\theta }_1,{\theta }_2)$ partitioned in such a way that there is a closed-form solution for $\pi ({\theta }_1|z,{\theta }_2)$ and $\pi ({\theta }_2|z,{\theta }_1)$ . For a given starting value for ${\theta }_1$ , say ${\theta }_1^0$ , we can set up an iterative procedure that allows us to draw a value for ${\theta }_2$ (denoted ${\theta }_2^1)$ given ${\theta }_1={\theta }_1^0$ . We then obtain a draw of ${\theta }_1$ (denoted ${\theta }_1^1)$ for the given ${\theta }_2={\theta }_2^1$ . Proceeding iteratively, we obtain a sequence of numerical draws for $\theta =({\theta }_1,{\theta }_2)$ . After the effect of initial conditions wears off, in many situations it can be shown that these draws are draws from the desired joint density $\pi (\theta |z)$ . Hence, after discarding a set of early draws, we can save the latter draws and use them to numerically characterize the posterior distribution for the parameters.

Example 5.3.1 (Least squares regression and the Gibbs sampler). Suppose $y_{t+1}$ is independently normally distributed conditional on $x_t$ and $\theta$ with mean ${\beta }^{\prime }{x}_t$ and variance $v^{-1}{I}_k,$ where $k=\dim (x_t)$ and $\theta =(\beta ,\upsilon )$ , both of which are unknown. Suppose we set the prior for the conditional mean $\beta \sim \mathrm{N}({\beta }_0,{\tau }^{2}I)$ independent of $y_{t+1}$ and the prior for the inverse of the variance, so v is a ${\chi }_{\omega }^2$ divided by a constant, $s_0,$ and the mean of the prior for the variance becomes $s_0/\omega$ . The hyperparameters are then ${\theta }_0=({\beta }_0,{\tau }^2,\omega ,s_0)$ . For this problem there is no simple analytical or closed-form solution for $\pi (\theta |z).$ , and numerical methods are required. It can be shown (Geweke, 2005, example 2.1.2) that

$$\beta |z,v,{\theta }_0\sim \mathrm{N}(\tilde{\beta }(v),W(v)),$$

$$s^2(\beta )v\sim {\chi }_{\tilde{\omega }}^2,$$

where $\begin{array}{r}W(v)=({\tau }^{-2}{I}_k+v{\sum }_{t=0}^{T-1}{x}_t{x}_t^{\prime }),\tilde{\beta }(v)=W^{-1}\left({\tau }^{-2}{\beta }_0+v{\sum }_{t=0}^{T-1}{x}_t{y}_{t+1}\right)\end{array}$ $\tilde{\omega }=\omega +T$ and $\begin{array}{r}{s}^2(\beta )=s_0^2+{\sum }_{t=0}^{T-1}(y_{t+1}-\beta x_t{)}^2\end{array}$ . Since it is straightforward to draw from the normal and ${\chi }^2$ distributions, we can easily compute a sequence of values from these conditional distributions.

Using draws of θ generated from $\pi (\theta |z)$ , it is straightforward to either directly construct the forecast when this is a known function of $\pi (\theta |z)$ , or alternatively construct draws of the forecast, $\begin{array}{r}{p}^P(y|z,\theta )\end{array}$ . In the first case the forecast $f(z)$ is a known function of features of the distribution $p(y|z,\theta )$ and of the distribution for $\theta .$ . For example, the linear regression case requires only the posterior means of θ given z which are easily computed. These can be used as the “coefficient estimates” in constructing the linear forecasting model.

The forecast density, ${\mathcal{P}}_y^P(y|z)$ , can be constructed from the draws of π $\left(\theta |z\right)$ and the known model, ${\varphi }_y(y|z,\theta )$ , drawing a value for the prediction for each of the draws of $\theta ,$ , while holding z constant. As noted in chapter 3, this approach can be used for any loss function, with different features of the estimated forecast distribution being relevant for different loss functions.

Suppose we are interested in the actual forecast rather than the predictive density, but there is no closed-form solution for the forecast as a function of this density. In this case we need to compute $f^{\ast }(z)=\arg$ min $r(\pi ,f(z))$ . Conditional on $z,$ this requires finding the rule that minimizes $\int L(f(z),y)p_Y^p(y|z)dy$ . Using the draws from $p_Y^P(y|z)$ , for any forecast rule, $f(z)$ , we can compute the average loss.

This follows because for any such forecast rule and any draw for y, we can compute a value for $L(f(z),y)$ with the average of these draws being the expected loss. Thus we can directly compute the expected loss associated with any forecasting rule. From among a set of such rules we can then choose the best rule as that with the lowest expected loss.

Zellner (1971) gives an early exposition of the Bayesian analysis of regression models including many of the models typically used in forecasting. More recent treatments and discussions of algorithms and numerical methods are provided in Koop (2003), Geweke (2005), and Karlsson (2013).