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.
79,328 characters · 21 sections · 82 citation commands
Bayesian Optimization of Hyperparameters from Noisy Marginal Likelihood Estimates
Keywords: acquisition strategy, Bayesian optimization, importance sampling, MCMC, steady-state BVAR, stochastic volatility, US macro.
The trend in econometrics is to use increasingly more flexible models that give a richer description of the economy, particularly for prediction purposes. As the model complexity increases, the estimation problems get more involved, and computationally costly MCMC methods are often used to sample from the posterior distribution.
Most models involve a relatively small set of hyperparameters that needs to be chosen by the user. For example, consider the steady-state BVAR model villani2009, which is widely used among practitioners and professional forecasters karlsson2013, and used in Section (ref) for illustration. The choice of the prior distribution in BVARs is often reduced to the selection of a small set of prior hyperparameters. Some of these hyperparameters can be specified subjectively by experts, for example, the steady-state is usually given a rather informative subjective prior. Other prior hyperparameters control the smoothness/shrinkage properties of the model and are less easy to specify subjectively.
giannone2015 proposed to treat these hard-to-specify prior hyperparameters as unknown parameters and explore the joint posterior of the hyperparameters, the VAR dynamics, and the shock covariance matrix. This is a statistically elegant approach which works well when the marginal likelihood is available in closed form and is easily evaluated. However, the marginal likelihood is rarely available in closed form. The BVARs with conjugate priors considered in carriero2012, and giannone2015 are an exception, but already the steady-state VAR needs MCMC methods to evaluate the marginal likelihood. It is of course always an option to sample the hyperparameters jointly with the other model parameters using Metropolis-Hastings (MH) or Hamiltonian Monte Carlo (HMC), but this likely leads to inefficient samplers since the parameter spaces are high-dimensional and the posterior of the hyperparameters are often quite complex, see e.g. the application in Section (ref).
Most practitioners also seem to prefer to fix the hyperparameters before estimating the other model parameters. carriero2012 propose a brute force optimization approach where the marginal likelihood is evaluated over a grid. This is computationally demanding, especially if a simulation based method has to be used for computing the marginal likelihood. Since the marginal likelihood in giannone2015 is available in closed form, one can readily optimize it using a standard gradient based optimizer with automatic differentiation, but this is again restricted to models with conjugate priors. The vast majority of applications instead use so-called conventional values for the hyperparameters, dating back to doan1984forecasting, which were found to be optimal on a specific historical dataset but are likely to be suboptimal for other datasets. Hence, there is a real need for a fast method for optimizing the marginal likelihood over a set of hyperparameters when every evaluation of the marginal likelihood is a noisy estimate from a computationally costly full MCMC run.
Bayesian optimization (BO) is an iterative optimization technique originating from machine learning. BO is particularly suitable for optimization of costly noisy functions in small to moderate dimensional parameter spaces (brochu2010 and snoek2012practical) and is therefore well suited for marginal likelihood optimization. The method treats the underlying objective function as an unknown object that can be inferred by Bayesian inference by evaluating the function at a finite number of points. A Gaussian process prior expresses Bayesian prior beliefs about the underlying function, often just containing the information that the function is believed to have a certain smoothness. Bayes theorem is then used to sequentially update the Gaussian process posterior after each new function evaluation. Bayesian optimization uses the most recently updated posterior of the function to decide where to optimally place the next function evaluation. This so-called acquisition strategy is a trade-off between: i) exploiting the available knowledge about the function to improve the current maxima and ii) exploring the function to reduce the posterior uncertainty about the objective function.
Our paper proposes a novel framework for Bayesian optimization when the user can control the precision and computational cost of each function evaluation. The framework is quite general, but we focus mainly on the situation when the noisy objective function is a marginal likelihood computed by MCMC. This is a very common situation in econometrics using, for example, the estimators in chib1995marginal, chib2001marginal, and geweke1999using. The precision of the marginal likelihood estimate at each evaluation point is then implicitly chosen by the user via the number of MCMC iterations. This makes it possible to use occasional cheap noisy evaluations of the marginal likelihood to quickly explore the marginal likelihood over hyperparameter space during the optimization. Our proposed acquisition strategy can be seen as jointly deciding where to place the new evaluation but also how much computational effort to spend in obtaining the estimate. We implement this strategy by a stopping rule for the MCMC sampling combined with an auxiliary prediction model for the computational effort at any new evaluation point; the auxiliary prediction model is learned during the course of the optimization.
We apply the method to the steady-state BVAR villani2009 and the time-varying parameter BVAR with stochastic volatility chan2018 and demonstrate that the new acquisition strategy finds the optimal hyperparameters faster than traditionally used acquisition functions. It is also substantially faster than a grid search and finds a better optimum.
The outline of the paper is a follows. Section (ref) introduces the problem of inferring hyperparameters from an estimated marginal likelihood. Section (ref) gives the necessary background on Gaussian processes and Bayesian optimization and introduces our new Bayesian optimization framework. Section (ref) illustrates and evaluates the proposed algorithm in a simulation study. Sections (ref) and (ref) assess the performance of the proposed algorithm in empirical examples, and the final section concludes.
To introduce the problem of using an estimated marginal likelihood for learning a model's hyperparameters we consider the selection of hyperparameters in the popular class of Bayesian vector autoregressive models (BVARs) as our running example.
Consider the standard BVAR model
where $\left\{ \boldsymbol{\varepsilon}_{t}\right\} _{t=1}^{T}$ are iid $N(\boldsymbol{0},\boldsymbol{\Sigma})$. A simplified version of the Minnesota prior (see e.g. karlsson2013) without cross-equation shrinkage is of the form
with
where $\text{vec}(\underline{\boldsymbol{\Pi}})$ and $\underline{\boldsymbol{\Omega}}_{\Pi}$ denotes the prior mean and covariance of the coefficient matrix, $\underline{\boldsymbol{S}}$ is the prior scale matrix with the prior degrees of freedom, $\underline{\nu}$. The diagonal elements of $\underline{\boldsymbol{\Omega}}_{\Pi}$ are given by
where $\lambda_{1}$ controls the overall shrinkage and $\lambda_{3}$ the lag-decay shrinkage set by the user, $s_{r}$ denotes the estimated standard deviation of variable $r$. The fact that we do not use the additional cross-equation shrinkage hyperparameter, $\lambda_{2}$, makes this prior conjugate to the VAR likelihood, a fact that will be important in the following. It has been common practice to use standard values that dates back to doan1984forecasting, but there has been a renewed interest to find values that are optimal for the given application (see e.g. banbura, carriero2012 and giannone2015). Two main approaches have been proposed. First, giannone2015 proposed to sample from the joint posterior distribution using the decomposition
where $\boldsymbol{\beta}=(\Pi,\Sigma)$ and $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ is the marginal posterior distribution of the hyperparameters. The algorithm samples from $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ using Metropolis-Hastings (MH) and then samples directly from $p(\boldsymbol{\beta}\vert\boldsymbol{\theta},\mathbf{y}_{1:T})$ for each $\boldsymbol{\theta}$ draw by drawing $\Pi$ and $\Sigma$ from the Normal-Inverse Wishart distribution. There are some limitations to using this approach. First, the $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ can be multimodal (see e.g. the application in Section (ref)) and it can be hard to find a good MH proposal density, making the sampling time-consuming. Second, practitioners tend to view hyperparameter selection as similar to model selection and want to determine a fixed value for $\boldsymbol{\theta}$ once and for all early in the model building process.
carriero2012 propose an exhaustive grid search to find the $\boldsymbol{\theta}$ that maximizes $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ and then uses that optimal $\boldsymbol{\theta}$ throughout the remaining analysis. The obvious drawback here is that a grid search is very costly, especially if we have non-conjugate priors and more than a couple of hyperparameters.
A problem with both the approach in giannone2015 and carriero2012 is that for most interesting models $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ is not available in closed form. Even the Minnesota prior with cross-equation shrinkage is no longer a conjugate prior and $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ is intractable. In fact, most Bayesian models used in practice have intractable $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$, including the steady-state BVAR villani2009 and the TVP-SV BVAR primiceri2005 used in Section (ref) and (ref).
When $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ is intractable, MCMC or other simulation based methods like Sequential Monte Carlo (SMC) are typically used to obtain a noisy estimate of $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$. Since
where $p(\mathbf{y}_{1:T}\vert\boldsymbol{\theta})=\int p(\mathbf{y}_{1:T}\vert\boldsymbol{\theta},\beta)p(\beta\vert\boldsymbol{\theta})d\beta$ is the marginal likelihood, this problem goes under the heading of (log) marginal likelihood estimation. We will therefore frame our method as maximizing the log marginal likelihood; maximization of $p(\boldsymbol{\theta}\vert\mathbf{y}_{1:T})$ is achieved by simply adding the log prior, $\log p(\boldsymbol{\theta})$, to the objective function.
chib1995marginal proposes an accurate way of computing a simulation-consistent estimate of the marginal likelihood when the posterior can be obtained via Gibbs sampling, which is the case for many econometric models. We refer to chib1995marginal for details about the marginal likelihood estimator and its approximate standard error.
Estimating the marginal likelihood for the TVP-SV BVAR is more challenging and we will adopt the method suggested by chan2018. The approach consists of four steps; 1) obtain a posterior sample via Gibbs sampling, 2) integrate out the time-varying VAR coefficients analytically, 3) integrate out the stochastic volatility using importance sampling to obtain the integrated likelihood, 4) integrate out the static parameters using another importance sampler. Since the algorithm makes use of two nested importance samplers, it is a special case of importance sampling squared ($IS^{2}$) tran2013importance; see Section (ref) for more details.
There are many alternative estimators that can be used in our approach, for example, the extension of Chib's estimator to Metropolis-Hastings sampling chib2001marginal, and estimators based on importance sampling geweke1999using or Sequential Monte Carlo doucet2001sequential.
All simulation-based estimators: i) give noisy evaluations of the marginal likelihood, ii) are time-consuming, and iii) have a precision that is controlled by the user in terms of the number of MCMC or importance sampling draws. The next section explains how traditional Bayesian optimization is well suited for points i) and ii), but lacks a mechanism for exploiting point iii). Taking the user controlled precision into account brings a new perspective to the problem and we propose a class of algorithms that handle all three points above.
Since Bayesian optimization is a relatively unknown method in econometrics, we give an introduction here to Gaussian processes and their use in Bayesian optimization.
A Gaussian process (GP) is a (possibly infinite) collection of random variables such that any subset is jointly distributed according to a multivariate normal distribution, see e.g. williams2006gaussian. This process, denoted by $f(\mathbf{x})\sim\mathcal{GP}(\mu(\mathbf{x}),k(\mathbf{x},\mathbf{x}'))$, can be seen as a probability distribution over functions $f:\mathcal{X}\to\mathbb{R}$ that is completely specified by its mean function, $\mu(\mathbf{x})\equiv\mathbb{E}f(\mathbf{x})$, and its covariance function, $\mathbb{C}(f(\mathbf{x}),f(\mathbf{x}'))\equiv k(\mathbf{x},\mathbf{x}')$, where $\mathbf{x}$ and $\mathbf{x}'$ are two arbitrary input values to $f(\cdot)$. Note that the covariance function specifies the covariance between any two function values, $f(\mathbf{x}_{1})$ and $f(\mathbf{x}_{2})$. A popular covariance function is the squared exponential (SE):
where $|\mathbf{x}-\mathbf{x}'|$ is the Euclidean distance between the two inputs; the covariance function is specified by its two kernel hyperparameters, the scale parameter $\sigma_{f}>0$ and the length scale $\ell>0$. The scale parameter $\sigma_{f}$ governs the variability of the function and the length scale determines how fast the correlation between two function values taper off with the distance $|\mathbf{x}-\mathbf{x}'|$, see Figure (ref). The fact that any finite sampling of function values $\lbrace f(\mathbf{x}_{n})\text{ for }\mathbf{x}_{n}\in\mathcal{X}\rbrace_{n=1}^{N}$ constitutes a multivariate normal distribution on $\mathbb{R}^{N}$ allows for the convenient conditioning and marginalization properties of the multivariate normal distribution. In particular, this makes it easy to compute the posterior distribution for the function $f$ at any input $\mathbf{x}_{\star}$.
An increasingly popular alternative to the squared exponential kernel is the Mat�rn kernel, see e.g. matern1960 and williams2006gaussian. The Mat�rn kernel has an additional hyperparameter, $\nu$>0, in addition to the length scale $\ell$ and scale $\sigma_{f}$, such that the process is $k$ times mean square differentiable if and only if $\nu>k$. Hence, $\nu$ controls the smoothness of the process and it can be shown that the Mat�rn kernel approaches the SE kernel as $\nu\rightarrow\infty$ (williams2006gaussian). Our approach is directly applicable for any valid kernel function, but we will use the popular Mat�rn $\nu=5/2$ kernel in our applications:
where $r=|\mathbf{x}-\mathbf{x}'|$. The Mat�rn 5/2 has two continuous derivatives which is often a requirement for Newton-type optimizers snoek2012practical. The hyperparameters, $\sigma_{f}$ and $\ell$, are found by maximizing the marginal likelihood, see e.g. williams2006gaussian.
Consider the nonlinear/nonparametric regression model with additive Gaussian errors
and the prior $f(\mathbf{x})\sim\mathcal{GP}(0,k(\mathbf{x},\mathbf{x}'))$. Given a dataset with $n$ observations, the posterior distribution of $f(\boldsymbol{x}_{\star})$ at a new input $\boldsymbol{x}_{\star}$ is williams2006gaussian
where $\boldsymbol{y}=(y_{1},\ldots,y_{n})^{\top}$, $k(\boldsymbol{x}_{\star})$ is the $n$-vector with covariances between $f$ at the test point $\boldsymbol{x}_{\star}$ and all other training inputs, $\boldsymbol{K}(\boldsymbol{X},\boldsymbol{X})$ is the $n\times n$ matrix with covariances among the function values at all $n$ training inputs in $\boldsymbol{X}$. When the errors are heteroscedastic with variance $\sigma_{i}^{2}$ for the $i$th observation, the same formulas apply with $\sigma^{2}I$ replaced by $\mathrm{diag}(\sigma_{1}^{2},\ldots,\sigma_{n}^{2})$.
Bayesian optimization (BO) is an iterative optimization method that selects new evaluation points using the posterior distribution of $f$ conditional on the previous function evaluations. More specifically, BO uses an acquisition function, $a(\boldsymbol{x})$, to select the next evaluation point (brochu2010 and snoek2012practical).
An intuitively sensible acquisition rule is to select a new evaluation point that maximizes the probability of obtaining a higher function value than the current maximum, i.e. the Probability of Improvement (PI):
where $f_{max}$ is the maximum value of the function obtained so far. The functions $m(\mathbf{x})$ and $s(\mathbf{x})$ are the posterior mean and standard deviation of the estimated Gaussian process for $f$ in the point $\mathbf{x}$, conditional on the available function evaluations (see ((ref))), and $\Phi$ denotes the cumulative standard normal distribution.
The Expected Improvement (EI) takes also the size of the improvement into consideration:
where $\phi$ denotes the density function of the standard normal distribution. The first part of ((ref)) is associated with the size of our predicted improvement and the second part is related to the uncertainty of our function in that area. Thus, EI incorporates the trade-off between high expected improvement (exploitation) and learning more about the underlying function (exploration). Optimization of the acquisition function, to find the next evaluation point, is typically a fairly easy task since it is noise-free and cheap to evaluate, however, it can have multiple optima so some care has to be taken. An easily implemented solution is to use a regular Newton-type algorithm initiated with different starting values over the hyperparameter surface. In this paper we use particle swarm optimization, a global optimization algorithm implemented in the Optim.jl package in Julia.
The exploitation-exploration trade-off is illustrated in Figure (ref), where the blue line shows the true objective function, the black line denotes the posterior mean of the GP, and the blue-shaded regions are 95% posterior probability bands for $f$. The black (small) dots are past function evaluations, and the violet (large) dot is the current evaluation. The first and third row in Figure (ref) show the true objective function and the learned GP at four different iterations of the algorithm while the second and fourth row show the EI acquisition function corresponding to the row immediately above. At Iteration 2 in the top-left corner, we see that the EI acquisition function (second row) indicates that there is a high expected improvement by moving to either the immediate left or right of the current evaluation. At Iteration 5 the EI strategy suggests three regions worth evaluating, where the two leftmost regions are candidates because of their high uncertainty. After seven iterations, the algorithm is close to the global maximum and will now continue a narrow search for the exact location of the global maximum.
Acquisition rules like PI or EI do not consider that different evaluation points can be more or less costly. To introduce the notation of cost into the acquisition strategy, snoek2012practical proposed Expected Improvement per second, $\text{EIS}(\mathbf{x})\equiv\text{EI}(\mathbf{x})/c(\mathbf{x})$, where $c:\mathcal{X}\to\mathbb{R}^{+}$ is a duration function that measures the evaluation time at input $\mathbf{x}$ in seconds. More generally, we can define $a(\mathbf{x})/c(\mathbf{x})$ as an effort-aware acquisition function. The duration function is typically unknown and snoek2012practical proposed to estimate it alongside $f$ using an additional Gaussian process for $\log c(\mathbf{x})$.
EIS assumes that the duration (or the cost) of function evaluations are unknown, but fixed for a given input $\mathbf{x}$; once we visit $\mathbf{x}$, the cost of the function estimate $\hat{f}(\mathbf{x})$ is given. However, the user can often choose the duration spent to obtain a certain precision in the estimate; for example by increasing the number of MCMC iterations when the marginal likelihood is estimated by MCMC. This novel perspective opens up for strategies that not only optimize for the next evaluation point, but also optimize over the computational resources, or equivalently, the precision of the estimate $\hat{f}(\mathbf{x})$. We formally extend BO by modeling the function evaluations with a heteroscedastic GP
where the noise variance $\sigma^{2}(\mathbf{x},G)$ is now an explicit function of the number of MCMC iterations, $G$, or some other duration measure. Hence the user can now choose both where to place the next evaluation and the effort spent in computing it by maximizing
with respect to both $\mathbf{x}$ and $G$, where $a(\mathbf{x})$ is a baseline acquisition function, for example EI.
A complication with maximization of $\tilde{a}(\mathbf{x},G)$ is that while we typically know that $\sigma(\mathbf{x},G)=O(1/\sqrt{G})$ in Monte Carlo or MCMC, the exact numerical standard error depends on the integrated autocorrelation time (IACT) of the MCMC chain. Note that the evaluation points can, for example, be hyperparameters in the prior, where different values can give rise to varying degrees of well-behaved posteriors, so we can not expect the IACT to be constant over the hyperparameter space, hence the explicit dependence on $\boldsymbol{x}$ in $\sigma^{2}(\mathbf{x},G)$. Rather than maximizing $\tilde{a}(\mathbf{x},G)$ with respect to both $\mathbf{x}$ and $G$ directly, we propose to implement the algorithm in an alternative way that achieves a similar effect. The approach includes stopping the evaluation early whenever the function evaluation turns out to be hopelessly low with a low probability of improvement over the current $f_{max}$.
For a given $\mathbf{x}$ we let $G$ increase, in batches of a fixed size, until
for some small value, $\alpha$, or until $G$ reaches a predetermined upper bound, $\bar{G}$; here $\hat{m}^{(g)}(\mathbf{x})\text{ and }s^{(g)}(\mathbf{x})$ denotes the posterior mean and standard deviation of the GP evaluated at $\mathbf{x}$ after $g$ MCMC iterations. Note that both the posterior mean $m(\mathbf{x})$ and standard deviation $s(\mathbf{x})$ are functions of the noise variance, which in turn is a function of $G$. The posterior distribution for $f(\mathbf{x})$ is hence continuously updated as $G$ grows until $1-\alpha$ of the posterior mass in the GP for $f(\mathbf{x})$ is concentrated below $f_{max}$, at which point the evaluation stops. The optimization is insensitive to the choice of $\alpha$, as long as it is a relatively small number. We now propose to maximize the following acquisition function based on early stopping
where $\hat{G}_{\alpha}(\mathbf{x})$ is a prediction of the number of MCMC draws needed at $\mathbf{x}$ before the evaluation stops, with the probability $\alpha$ as the threshold for stopping. We emphasize that early stopping is here used in a subtle way, not only as a simple rule to short-circuit useless computations, but also in the planning of future computations; the mere possibility of early stopping can make the algorithm try an $\mathbf{x}$ which does not have the highest $a(\mathbf{x})$, but which is expected to be cheap and is therefore worth a try. This effect that comes via $\sigma^{2}(G)$ is not present in the EIS of snoek2012practical where the cost is fixed and is not influenced by the probability model on $f$.
Although one can use any model to predict $G$, we will here fit a GP regression model to the logarithm of the number of MCMC draws, $\log G_{j}$ for $j=1,...,J$ in the $J$ previous evaluations
where $\mathbf{z}_{j}$ is a vector of covariates. The hyperparameters, $\mathbf{x}_{1:J}$, themselves may be used as predictors of $\hat{G}(\mathbf{x})$, but also $D(\mathbf{x})=\hat{m}(\mathbf{x})-f_{max}$ and $s(\mathbf{x})$ are likely to have predictive power for $G$, as well as $u(\mathbf{x})=(\hat{m}(\mathbf{x})-f_{max})/s(\mathbf{x})$. We will use $\mathbf{z}_{j}=\left(\mathbf{x}_{j},D^{(j)}(\mathbf{x}_{j}),s^{(j)}(\mathbf{x}_{j}),u^{(j)}(\mathbf{x}_{j})\right)$ in our applications, where the superscript over $j$ denotes the BO iteration. The prediction for $G$ is taken to be $\hat{G}=\text{exp}\left(m_{G}\left(\mathbf{z}\right)\right)$, which corresponds to the median of the log-normal posterior for $G$.
We will use the term Bayesian Optimization with Optimized Precision (BOOP) for BO methods that optimize $\tilde{a}_{\alpha}(\mathbf{x})$ in ((ref)), and more specifically BOOP-EI when EI is used as the baseline acquisition function, $a(\mathbf{x})$. The whole procedure is described in Algorithm (ref).
Note that ((ref)) assumes that $\hat{f}(\mathbf{x})$ is an unbiased estimator at any $\mathbf{x}$. This can be ensured by using enough MCMC/importance sampling draws in the first marginal likelihood evaluation of BOOP. We performed a small simulation exercise that shows that the Chib estimator is approximately unbiased after a small number of iterations for the medium-sized VAR model in (ref). As expected, we had to use more initial draws in the large-scale VAR in Section (ref) and the time-varying parameter VAR with stochastic volatility in Section (ref) to drive down the bias. See also Section (ref) for some ideas on how to extend BOOP to estimators where the bias is still sizeable in large MCMC samples.
Figure (ref) illustrates the early stopping part of BOOP in a toy example. The first row illustrates the first BOOP iteration and the columns show increasingly larger MCMC sample sizes ($G$). We can see that the 95% posterior interval after $G=10$ MCMC draws at the current $\mathbf{x}$ includes $f_{max}^{(1)}$ (dotted orange line), the highest posterior mean of the function values observed so far; it therefore worthwhile to increase the number of simulations for this $\mathbf{x}$. Moving one graph to the right we see that after $G=20$ simulations the 95% posterior interval still includes $f_{max}^{(1)}$, and we move one more graph to the right for $G=50$. Here we conclude that the sampled point is almost certainly not an improvement and we move on to a new evaluation point. The new evaluation point is found by maximizing the BOOP-EI acquisition function in ((ref)) with updated effort prediction function $\hat{G}(\boldsymbol{z})$ in Equation (ref) and is depicted by the violet dot in leftmost graph in the second row of Figure (ref). Following the progress in the second row, we see that it takes only $G=20$ samples to conclude that the function value is almost certainly lower than the current maximum of the posterior mean at the second BO iteration. Finally, in the third row, we can see that the point is sampled with high variance at the beginning, but as we increase $G$ it becomes clear that this $\mathbf{x}$ is indeed an improvement.
The code used in this paper is written in the Julia computing language bezanson2017julia, making use of the GaussianProcesses.jl package for estimation of the Gaussian processes and the Optim.jl package for the optimization of the acquisition functions.
We will here assess the performance of the proposed BOOP-EI algorithm in a simulation experiment for the optimization of a non-linear function in a single variable. The simple one-dimensional setup is chosen for presentational purposes, and more challenging higher-dimensional settings are likely to show even larger advantages of our method compared to regular Bayesian optimization.
The function to be maximized is
where $N(x|\mu,\sigma^{2})$ denotes the density function of a $N(\mu,\sigma^{2})$ variable; the function is plotted in Figure (ref) (left). We further assume that $f(x)$ can be estimated at any $x$ from $G$ noisy evaluations by a simple Monte Carlo average
where $g(x)=0.1+0.15N(x|1,0.5^{2})+0.15N(x|2.5,0.25^{2})+0.5N(x|5,0.75^{2})$ is a heteroscedastic standard deviation function that mimics the real world case where the variability of the marginal likelihood estimate varies over the space of hyperparameters; $g(x)$ is plotted in Figure (ref) (middle). Figure (ref) (right) illustrates the noisy function evaluations (gray points) and the effect of Monte Carlo averaging $G=3$ evaluations for a given $x$ (blue points). We assume for simplicity here that once the algorithm decides to visit an $x$ it will get access to a noise-free evaluation of the standard error of the estimator.
Before we conduct a small Monte Carlo experiment it is illustrative to look at the results from a single run of the EI and BOOP-EI algorithms. Figure (ref) highlights the difference between the algorithms by showing the GP posterior after 25 Bayesian optimization runs. The EI algorithm has clearly wasted computations to needlessly drive down the uncertainty at (useless) low function values, while BOOP-EI tolerates larger uncertainty in such function locations.
We now compare the performance of a standard EI approach using a heteroscedastic GP with our BOOP-EI approach in a small simulation study. The methods will be judged by their ability to find the optimum using as few evaluations as possible. The performance will be evaluated by a Monte Carlo study where we simulate $1000$ replications from each model under each simulation scenario. We investigate Bayesian optimization with EI using $100$ and $500$ samples in each iteration, BOOP-EI is allowed to stop the sampling at any time before the number of samples for the EI is reached. In each simulation we set an upper bound of $120$ Bayesian optimization iterations.
We can see from Figure (ref) that BOOP finds the maximum by using fewer total number of samples in both scenarios and that the difference increases with the number of samples used when forming the estimator. We can also see from the median $f_{max}$ that both algorithms can get stuck for a while at the second greatest local maximum (which is approximately $.05$ lower than the global maximum). However, BOOP gets out of the local optimum faster since it has the option to try cheap noisy evaluations and will therefore explore other parts of the function earlier than basic BO-EI. This effect seems to increase as we allow for a higher number of samples since it lowers the relative price of cheaper evaluations.
In this section, we use BOOP-EI to estimate the prior hyperparameters of the steady-state BVAR of villani2009. giannone2015 show that finding the right values for the hyperparameters in BVARs can significantly improve forecasting performance. Moreover, banbura show that different degree of shrinkage (controlled by the hyperparameters) is necessary under different model specifications.
The steady-state BVAR model of villani2009 is given by:
where $E[\mathbf{y}_{t}]=\Psi\mathbf{x}_{t}$. In particular, if we assume that $\mathbf{x}_{t}=1$ for all $t,$ then $\Psi$ has the interpretation as the overall mean of the process. We take the prior distribution to be:
where $\underline{\boldsymbol{\theta}}_{\Psi}\text{ and }\underline{\Omega}_{\Psi}$ are the mean and covariance matrix for the steady states. The prior covariance matrix for the dynamics, $\underline{\Omega}_{\Pi}$, is constructed using
where $\underline{\omega}_{ii}$ is the diagonal elements of $\underline{\Omega}_{\Pi}$. We also assume prior independence, following villani2009. The hyperparameters that we optimize over are: the overall-shrinkage parameter $\theta_{1}$, the cross-lag shrinkage $\theta_{2}$, and the lag-decay parameter $\theta_{3}$.
The posterior distribution of the steady-state BVAR model parameters can be sampled with a simple Gibbs sampling scheme villani2009. The marginal likelihood, together with its empirical standard error, can be estimated by the method in chib1995marginal.
Table (ref) describes the data used in our applications which are also used in giannone2015. It contains 23 macroeconomic variables for which two subsets are selected to represent a medium-sized model with 7 variables and a large model that contains 22 of the variables (real investment is excluded). Before the analysis, the consumer price index and the five-year bond rate were transformed from monthly to quarterly frequency. All series are transformed such that they become stationary according to the augmented Dickey-Fuller test. This is necessary for the data to be consistent with the prior assumption of a steady-state. The number of lags is chosen according to the HQ-criteria, hannan1979determination and quinnMulti. This resulted in $p=2$ lags for the medium-sized model which we also use for the large model.
We set the prior mean of the coefficient matrix, $\Pi$, to values that reflect some persistence on the first lag, but also that all the time series are stationary; e.g. the prior mean on the first lag of the FED interest rate and the GDP-deflator is set to 0.6, while others are set to zero in the medium-sized model. Lags longer than 1 and cross-lags all have zero prior means. The priors for the steady-states are set informative to the values listed in Table (ref), these values follow suggestions from the literature for most variables, see e.g. louzis2019 and osterholm2012limited. There were a few variables where we could not find theoretical values for either the mean or the standard deviation, in those cases, we set them close to their empirical counterparts.
We consider three competing optimization strategies: (I) an exhaustive grid-search, (II) Bayesian optimization with the EI acquisition function (BO-EI), and (III) our BOOP-EI algorithm. In each approach, we use the restrictions $\theta_{1}\in(0,5),\theta_{2}\in(0,1)\text{, and }\theta_{3}\in(0,5)$. In the grid-search, $\theta_{1}\text{ and }\theta_{2}$ move in steps of $0.05$ and $\theta_{3}$ moves in steps of $0.1$, yielding in total 100000 marginal likelihood evaluations. For the Bayesian optimization algorithm, we set the number of evaluations to 250, and we use three random draws as initial values for the GPs.
For strategies (I) and (II) we use a total of 10000 Gibbs iterations with 1000 as a burn-in sample in each model evaluation. For (III) we first draw 1100 Gibbs samples where we discard the first 1000 as burn-in and use the rest to calculate the probability of improvement PI, to ensure that the estimated marginal likelihood will be approximately unbiased; Figure (ref) shows that Chib's estimate is unbiased already after a few hundred samples. If $\text{PI}<\alpha$ we stop early and move on to the next BO iteration, otherwise we generate a new batch (of size 100) of Gibbs samples and again check the PI criteria. The total number of Gibbs iterations will therefore vary between 1100 and 10 000 in each of the 250 BOOP-iterations for the medium-sized model. Note that Chib's estimator uses an estimate of the parameter for the so-called reduced Gibbs sampler run. This point estimate should preferably have a high posterior density for efficiency reasons, see chib1995marginal. The medium-sized model uses only 100 posterior samples to obtain high-density parameters for calculating Chib's log marginal likelihood, which is enough in our set-up. For the large model, we use $5000$ burn-in samples and $500$ simulations in the first batch, which gives between $5500$ and $10000$ MCMC iterations per evaluation point. The $5000$ burn-in is likely to be excessive but is used to be conservative; a small number would make BOOP even faster in comparison to regular BO. The application is robust to the choice of $\alpha$, as long as it is a reasonably low number, in this study we use $\alpha=0.001$.
For comparison, we will also use the standard values of the hyperparameters used in e.g. the BEAR-toolbox, bear, $\theta_{1}=0.1,\theta_{2}=0.5,\text{ and }\theta_{3}=1$, as a benchmark. The methods are compared with respect to i) the obtained marginal likelihood, and ii) how much computational resources were spent in the optimization.
Table (ref) summarizes the results from ten runs of the algorithms for the medium size BVAR model. We see that all three optimization strategies find hyperparameters that yield substantially higher log marginal likelihood than the standard values. We can also see that both Bayesian optimization methods yield as good hyperparameters as the grid search at only a small fraction of the computational cost. It is also clear from Table (ref) that a substantial amount of computations associated with the MCMC are saved when using BOOP. It is interesting to note that the values for $\theta_{1}\text{ and }\theta_{2}$ are similar for all three optimization approaches but that $\theta_{3}$ differs to some extent. This is due to the flatness of the log marginal likelihood in that area.
The left graph of Figure (ref) shows that BOOP-EI finds higher values of the log marginal likelihood using much fewer MCMC iterations than plain BO with EI acquisitions. From Table (ref) we can see that BOOP-EI uses, on average, less than a fifth of the MCMC iterations compared to BO-EI for a full run. Interesting to note is that BO-EI leads to (on average) a higher number of improvements on the way to the maximum; while BOOP-EI gives fewer improvements but of larger magnitude; the strategy of cheaply exploring new territories before locally optimizing the function pays off. The graph to the right in Figure (ref) shows that for this application BO-EI is quicker in terms of CPU time to reach fairly high values for the log marginal likelihood. We see at least two explanations for this: first, BOOP-EI tries to explore more unknown territories since they are presumed to be cheap while BO-EI more greedily focuses on local optimization. Second, and more importantly, the overhead cost associated with the BOOP-EI acquisition is relatively large in this medium-sized application where the cost of evaluating the marginal likelihood itself is not excessive. The fact that BOOP can heavily reduce the number of MCMC draws while still giving similar CPU computational time suggest that it is most useful in cases where each log marginal likelihood evaluation is expensive; this will be demonstrated in the more computationally demanding models in Sections (ref) and (ref).
Figure (ref) displays the log marginal likelihood surfaces over the grid of $(\theta_{1},\theta_{2})$-values used in the grid search. Each sub-graph is for a fixed value of $\theta_{3}\in\{0.76,1,2\}$. The red dot indicates the predicted maximum log marginal likelihood for the given $\theta_{3},$ and the black dot in the middle sub-figure indicates the standard values. We can see that the standard values are located outside the high-density region, relatively far away from the maximum. A comparison of Figures (ref) and (ref) shows that the GP's predicted log marginal likelihood surface is quite accurate already after merely 250 evaluations; this is quite impressive considering that Bayesian optimization tries to find the maximum in the fastest way, and does not aim to have high precision in low-density regions.
We also optimize the parameters of the more challenging large BVAR model containing the 22 different time series, using 250 iterations for both BO-EI and BOOP-EI. A complete grid search is too costly here, so we instead compare with parameters obtained from BOOP in the medium-sized BVAR in Section (ref), which is a realistic strategy in practical work.
Table (ref) shows that our method, again, finds optimal hyperparameters with dramatically larger log ML than standard values, and also substantially better values than those that are optimal for the medium-scale BVAR. Finally, note that the hyperparameters selected by BOOP-EI in the large-scale BVAR are quite different from those in the medium-scale model. The optimal $\theta_{1}$ applies less baseline shrinkage than before, but the lag decay ($\theta_{3}$) is higher, and in particular, the cross-lag shrinkage, $\theta_{2}$, is much closer to zero, implying much harder shrinkage towards univariate AR-processes. This latter result strongly suggests that the computationally attractive conjugate prior structure is a highly sub-optimal solution since such a prior requires that $\theta_{2}=1$. We can see that for this more computationally demanding model BOOP-EI is much faster and finish on average in a third of the time than the regular BO-EI strategy. Figure (ref) show the predicted log marginal likelihood surface obtained from the last GP in a BOOP run. The rightmost graph conditions on $\theta_{3}=1.51$, which is optimal for BOOP-EI, so this graph therefore has the GP with the highest accuracy.
The time-varying parameter BVAR with stochastic volatility (TVP-SV BVAR) in chan2018 is given by
where $\boldsymbol{\mu}_{t}$ is a vector of time varying intercepts, $\boldsymbol{B}_{1t},\dots,\boldsymbol{B}_{pt}$ are $n\times n$ matrices of VAR coefficients, $\boldsymbol{B}_{0t}$ is a $n\times n$ lower triangular matrix with ones on the main diagonal. The evolution of $\boldsymbol{\Sigma}_{t}=\mathrm{diag}\left(\exp(h_{1t}),\dots,\exp(h_{nt})\right)$ is modelled by the vector of log volatilities, $\mathbf{h}_{t}=(h_{1t},\dots,h_{nt})^{\top}$, evolving as a random walk
where $\boldsymbol{\Sigma}_{h}=\mathrm{diag}(\sigma_{h_{1}}^{2},\dots,\sigma_{h_{n}}^{2})$ and the starting values in $\mathbf{h}_{0}$ are parameters to be estimated. Following chan2018, we collect all parameters of $\boldsymbol{\mu}_{t}$ and the $\mathbf{B}_{it}$ matrices in a $k_{\gamma}$-dimensional vector $\boldsymbol{\gamma}_{t}$ and write the model in state space form as
where $\boldsymbol{X}_{t}$ contain both current and lagged values of $\mathbf{y}$, $\boldsymbol{\Sigma}_{\gamma}=\mathrm{diag}(\sigma_{\gamma_{1}}^{2},\dots,\sigma_{\gamma_{k\gamma}}^{2})$ and the initial values for the state variables follow $\boldsymbol{\boldsymbol{\gamma}}_{0}\sim N(\mathbf{a}_{\boldsymbol{\gamma}},\mathbf{V}_{\boldsymbol{\gamma}})$ and $\mathbf{h}_{0}\sim N(\mathbf{a}_{\mathbf{h}},\mathbf{V}_{\mathbf{h}})$.
For comparability we choose to use the same prior setup as in chan2018 where the variances of the state innovations follow independent inverse-gamma distributions: $\sigma_{\gamma_{i}}^{2}\sim IG(\nu_{\gamma_{0}},S_{\gamma_{0}})$ if $\gamma_{i}$ is an intercept, $\sigma_{\gamma_{i}}^{2}\sim IG(\nu_{\gamma_{1}},S_{\gamma_{1}})$ if $\gamma_{i}$ is a VAR coefficient and $\sigma_{h}^{2}\sim IG(\nu_{h},S_{h})$ for the innovations to the log variances. Following chan2018 we set $\mathbf{a}_{\boldsymbol{\gamma}}=\boldsymbol{0}$, $\mathbf{V}_{\boldsymbol{\gamma}}=10\cdot I_{k_{\gamma}}$, $\mathbf{a}_{\mathbf{h}}=0$, $\mathbf{V}_{\mathbf{h}}=10\cdot I_{n}$, and $\nu_{\gamma_{0}}=\nu_{\gamma_{1}}=\nu_{h}=5$, but we use Bayesian optimization to find the optimal values for the three key hyperparameters $S_{\gamma_{0}}$, $S_{\gamma_{1}}$ and $S_{h}$ which controls the degree of time-variation in the states. We collect the three optimized hyperparameter in the vector $\boldsymbol{\theta}=(\theta_{1},\theta_{2},\theta_{3})^{\top}$, where $\theta_{1}=S_{\gamma_{0}}$, $\theta_{2}=S_{\gamma_{1}}$ and $\theta_{3}=S_{h}$. We optimize over the domain $\{\boldsymbol{\theta}:0\leq\theta_{1}\leq5,0\leq\theta_{2}\leq1,0\leq\theta_{3}\leq5\}$, which allows for all cases from no time variation in any parameter to high time variation in all the model parameters.
To estimate the marginal likelihood, chan2018 first obtain posterior draws of $\boldsymbol{\gamma},\boldsymbol{h},\Sigma_{\boldsymbol{\gamma}},\Sigma_{\mathbf{h}},\boldsymbol{\gamma}_{0},\mathbf{h}_{0}$ using Gibbs sampling, which are then used to design efficient importance sampling proposals. The marginal likelihood is
where $\boldsymbol{\psi}$ collects all the static parameters in $\Sigma_{\boldsymbol{\gamma}},\Sigma_{\mathbf{h}},\boldsymbol{\gamma}_{0},\mathbf{h}_{0}$. The inner integral w.r.t. $\boldsymbol{\gamma}$ can be solved analytically and afterwards, we can integrate out $\mathbf{h}$ using importance sampling to obtain an estimate of the integrated likelihood
The last step is to integrate out the fixed parameters from the integrated likelihood
which is done using another importance sampler. The two nested importance samplers put the algorithm in the framework of importance sampling squared ($IS^{2}$, tran2013importance). The Chan-Eisenstat algorithm is elegantly designed, but necessarily computationally expensive with a single estimate of the marginal likelihood taking 205 minutes in MATLAB on a standard desktop chan2018. We call their MATLAB code from Julia using the MATLAB.jl package, illustrating that BOOP can plug in any marginal likelihood estimator. However, we found that the standard errors in chan2018 can be more robustly estimated using the bootstrap and we have done so here. The cost of bootstrapping the standard errors only has to be taken once for every Bayesian optimization iteration, and this cost is negligible compared to the computation of the log marginal likelihood estimate.
We use quarterly data for the GDP-deflator, real GDP, and the short-term interest rate in the USA from 1954Q3 to 2014Q4 from chan2018 for comparability. In addition, we make use of their Matlab code (with minor adjustments) for computing the marginal likelihood. This shows another strength with the BOOP approach, that it works on top of existing code.
We fix the number of Gibbs sampling iterations and the burn-in period to $20000$ and $5000$ respectively for both BO-EI and BOOP-EI in all evaluation points. This simplifies the implementation and does not make a practical difference since the main part of the computational cost is spent on the log marginal likelihood estimate from importance sampling. For BO-EI we use $5000$ log marginal likelihood evaluations in each new evaluation point, while BOOP-EI starts with $1000$ importance sampling draw and then takes batches of size $100$ until a maximum of $5000$ samples has been reached. The initial 1000 draws were enough to make the estimator approximately unbiased.
Table (ref) shows the optimized log marginal likelihood from three independent runs of BO-EI and BOOP-EI; the hyperparameter values used in chan2018 are shown for reference. As expected both BO and BOOP find better hyperparameters than the ones in chan2018; this is particularly true for BOOP which gives an increase in the marginal likelihood of more than 10 units on the log scale on average. Interestingly, both BO and BOOP suggest that the stochastic volatilities should be allowed to move more freely than in chan2018, but that there should be less time variation in the intercepts and VAR coefficients. This points in the same direction as the results in chan2018 who find that shutting down the time variation in the intercept and VAR coefficients actually increases the marginal likelihood. Our results indicate that when carefully selecting the shrinkage parameters by optimization, the VAR-dynamics should in fact be allowed to evolve over time, but at a slower pace.
Table (ref) shows a great deal of variability between runs, in particular for BO. Figure (ref) shows that this is probably because the hyperparameter surface is substantially more complicated and multimodal than for the steady-state BVAR. We can also see that the log marginal likelihood is relatively insensitive to changes in $\theta_{3}$ around the mode while it is very sensitive to changes in $\theta_{2}$.
We propose a new Bayesian optimization method for finding optimal hyperparameters in econometric models. The method can be used to optimize any noisy function where the precision is under the control of the user. We focus on the common situation of maximizing a marginal likelihood evaluated by MCMC or importance sampling, where the precision is determined by the number of MCMC or importance sampling draws. The ability to choose the precision makes it possible for the algorithm to take occasional cheap and noisy evaluations to explore the marginal likelihood surface, thereby finding the optimum faster.
We assess the performance of the new algorithm by optimizing the prior hyperparameters in the extensively used BVAR with stochastic volatility and time-varying parameters and the steady-state BVAR model in both a medium-sized and a large-scale VAR. The method is shown to be practical and competitive to other approaches in that it finds the optimum using a substantially smaller computational budget, and has the potential of being part of the standard toolkit for BVARs. We have focused on optimizing the marginal likelihood, but the method is directly applicable to other score functions, e.g. the popular log predictive score (geweke2007smoothly and villani2012generalized).
Our approach builds on the assumption that the noisy estimates of the log marginal likelihoods are approximately unbiased, which we verify is a reasonable assumption in the three applications if the first BOOP evaluation is based on a marginal likelihood estimator from enough MCMC draws. The unbiasedness of the log marginal likelihood will, however, depend on the combination of MCMC sampler and marginal likelihood estimator, see adolfson2007bayesian for some evidence from Dynamic Stochastic General Equilibrium (DSGE) models an2007bayesian. For example, the simulations in adolfson2007bayesian suggest that sampling with the independence Metropolis-Hastings combined with the chib2001marginal estimator is nearly unbiased, whereas sampling with the random walk Metropolis algorithm combined with the modified harmonic estimator geweke1999using can be severely biased, unless the posterior sample is extremely large. It would therefore be interesting to extend the method to cases with biased evaluations where the marginal likelihood estimates are persistent and only slowly approaching the true marginal likelihood. Since the marginal likelihood trajectory over MCMC iterations is rather smooth adolfson2007bayesian one can try to predict its evolution and then correct the bias in the marginal likelihood estimates.