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.
91,711 characters · 16 sections · 19 citation commands
-0.5 in Sequential Monte Carlo With Model Tempering
{ {\em Key words:} Bayesian Computations, Dynamic Stochastic General Equilibrium Models, Sequential Monte Carlo, Stochastic Volatility, Vector Autoregressions.}
\thispagestyle{empty}
\setcounter{page}{1}
Modern macroeconometrics often relies on time series models for which it is time-consuming to evaluate the likelihood function, either because it takes a long time to solve the underlying structural model, or the likelihood evaluation requires to integrate out latent state variables. In this paper we demonstrate how Bayesian computations for such models can be accelerated by reweighting and mutating posterior draws from an approximating model that allows for fast likelihood evaluations. We show that a sequential Monte Carlo (SMC) algorithm that starts out with draws from the posterior distribution of an approximating model instead of the prior distribution of the target model can drastically speed up the posterior computations.
SMC methods have been traditionally used to solve nonlinear filtering problems, an example being the bootstrap particle filter of GordonSalmondSmith1993. Subsequently, Chopin2002 showed how to adapt particle filtering techniques to conduct posterior inference for a static parameter vector. The first paper that applied SMC techniques to posterior inference for the parameters of a (small-scale) DSGE model was Creal2007. Subsequent work by herbst2014sequential,HerbstSchorfheide2016 fine-tuned the algorithm so that it could be used for the estimation of medium- and large-scale models. DurhamGeweke2011 show how to parallelize a flexible and self-tuning SMC algorithm for the estimation of time series models on graphical processing units (GPU).
In general, SMC algorithms approximate a target posterior distribution by creating intermediate approximations to a sequence of bridge distributions, indexed in this paper by $n$. At each stage, the current bridge distribution is represented by a swarm of so-called particles. Each particle is composed of a value and a weight. Weighted averages of the particle values converge to expectations under the stage-$n$ distribution. The transition from stage $n-1$ to $n$ involves changing the particle weights and values (mutation) so that the swarm adapts to the new distribution. Typically, these bridge distributions are constructed by either using the full-sample likelihood ({\em likelihood tempering}, LT)---generated by raising this likelihood function to the power of $\phi_n$, where $\phi_n$ increases from zero to one---or by sequentially adding observations to the likelihood function ({\em data tempering}, DT).
As initially suggested in CaiEtAl2020 but not explored any further, in this paper we document the runtime reductions achievable by a {\em model tempering approach} that takes a geometric average with weights $\phi_n$ and $1-\phi_n$ of the likelihood functions associated with the target model, denoted by $M_1$, and an approximating model $M_0$.\footnote{AcharyaEtAl2021 discuss model tempering as a strategy to estimate a HANK model.} Building on earlier work in the statistics literature, e.g., JasraStephensDoucetTsagaris2011, and work in the DSGE model literature, e.g., herbst2019tempered and Cai {\em et al.} (2021), we choose the tempering schedule defined through the $\phi_n$ sequence adaptively. Our adaptive schedules are calibrated by a single tuning parameter that controls the desired variance of the particle weights. The smaller the discrepancy between the posterior distribution of the approximating and the original model, the fewer bridge distributions are being used, and the faster the posterior analysis.
In general, model tempering is an attractive computational strategy for applications in which the likelihood evaluation for the target model is computationally costly and there is an approximating model for which the likelihood evaluation is fast and generates a posterior that is not too different from the posterior of the target model. We envision the approximating model to be a simplified version of the target model for which posterior computations are also implemented via SMC, in this case with likelihood tempering.\footnote{Even in the absence of a model tempering strategy, estimating approximating models is desirable as part of the modeling and code debugging that ultimately leads to the target model.} The $M_0$ likelihood tempering can be terminated before the weight on the likelihood function has reached the value one. We denote the terminal weight on the $M_0$ likelihood by $\psi_* \in (0,1]$. The early termination will lead to a more diffuse $M_0$ posterior, draws from which might be more easily mutable into draws from the $M_1$ posterior in the subsequent model tempering steps. This feature introduces additional flexibility into the model tempering algorithm.
We provide a formula for the runtime reduction achievable by model tempering that depends on the number of stages as a function of $\psi_*$ used for the $M_0$ and $M_1$ SMC runs, respectively, and the relative time it takes to evaluate the likelihood functions of the two models, denoted by the ratio $\tau_0/\tau_1$. Note that the user can evaluate $\tau_0/\tau_1$ before running the entire algorithm. We show that the runtime reduction profile is convergent as $\tau_0/\tau_1 \longrightarrow 0$. In the limit, the runtime reduction is determined just by the number of $M_0$ and $M_1$ SMC stages, which in turn depends on the alignment of the $\psi_*$-tempered $M_0$ posterior and the target $M_1$ posterior, relative to the alignment of the prior and the $M_1$ posterior. To assess the potential gains of model tempering {\em ex ante}, we recommend that the researcher computes the variance of the importance sampling weights, that would be needed to reweight the draws from the $\psi_*$-tempered $M_0$ posterior to approximate the target $M_1$ posterior, for various choices of $\psi_*$. If there is a $\psi_*$ for which this variance is small relative to the number of SMC particles, then the gains from model tempering are potentially large.
We consider three numerical illustrations of model tempering. In the first illustration, both target and approximating densities are univariate Normal. We illustrate how the distance between the densities affects the number of stages (and computational time) required to convert draws from the approximating density into draws from the target density. In the second illustration we consider the estimation of a vector autoregression (VAR) with stochastic volatility (SV), using a homoskedastic VAR as approximating model. In our illustration, model tempering is able to reduce the computational time by 79%. At last we consider the estimation of a dynamic stochastic general equilibrium (DSGE) model. We take $M_1$ as a version of the model that is solved with a second-order perturbation around the steady state and for which the likelihood function is evaluated with a bootstrap particle filter (BSPF). The approximating model is a log-linearized version for which the likelihood function can be evaluated quickly using the Kalman filter. In our numerical example, model tempering can reduce the runtime of our JULIA code from 655 to 80 minutes.
The remainder of this paper is organized as follows. Section (ref) describes the proposed model tempering SMC algorithm. Section (ref) considers the simple example based on univariate Gaussian posterior distributions. In Section (ref) we use model tempering to estimate the VAR with SV. In Section (ref) we implement our algorithm for the nonlinear DSGE model. Section (ref) concludes. An Online Appendix contains supplemental information on the methodology and further details and results for the numerical illustrations.
VARs, DSGE models, and other time series models are often estimated using Bayesian inference for several reasons. First, the Bayesian framework provides a powerful toolkit to handle the presence of latent variables in state-space models. Second, uncertainty about parameters, shocks, and unobserved state variables is treated identically which makes it conceptually straightforward to form predictive distributions that reflect all sources of uncertainty. Third, prior distributions can be used to regularize the estimation of high-dimensional models (e.g., VARs) or to incorporate additional information not contained in the estimation sample (DSGE model estimation).
Bayesian inference combines a prior distribution $p(\theta)$ with a likelihood function $p(Y|\theta)$ to form a posterior distribution $p(\theta|Y)$, which is given by
where $Y=Y_{1:T}=\{y_1, y_2, ..., y_T\}$ and the normalization constant $p(Y)$ is called the marginal data density (MDD). In most applications, the posterior distribution $p(\theta|Y)$ does not belong to a family of distributions for which moments and percentiles can be easily calculated or draws can be obtained by direct sampling. In this paper we use an SMC algorithm to sample from the posterior distribution $p(\theta|Y)$. The algorithm combines insights from importance sampling and Markov chain Monte Carlo (MCMC) techniques. Two of its key advantages are that it is able to provide accurate approximations of non-regular posterior distributions and that it can be easily parallelized, unlike MCMC algorithms. In Section (ref) we describe a generic SMC algorithm to sample from the posterior distribution of $\theta$. The section draws heavily from the more detailed exposition in herbst2014sequential,HerbstSchorfheide2016. Model tempering, which is the focus of our paper, is introduced in Section (ref) and implementation details are discussed in Section (ref). In Section (ref) we assess potential runtime reductions.
In order to draw from $\pi(\theta)$, the SMC algorithm uses a sequence of bridge posterior distributions $\{\pi_n(\theta)\}_{n=0}^{N_\phi}$, illustrated in Figure (ref), where the last one in the sequence equals the posterior distribution -- $\pi_{N_\phi}(\theta)=\pi(\theta)$ -- and where each $\pi_{n-1}(\theta)$ is used as the proposal density for $\pi_n(\theta)$. The bridge posteriors are constructed from stage $n$ likelihood functions $p_n(Y|\theta)$ and defined as
Each density $\pi_n(\theta)$ is represented by a particle approximation $\{\theta^i_n,W^i_n\}_{i=1}^N$. Thus, at stage $n$ the algorithm propagates the particles $\{\theta^i_{n-1},W^i_{n-1}\}_{i=1}^N$ so that they come to represent the target density $\pi_n(\theta)$. Formally, the algorithm proceeds in the following steps:
The correction step is a classic importance sampling step, in which the particle weights are updated to reflect the stage $n$ distribution $\pi_n(\theta)$. The selection step is optional. On the one hand, resampling adds noise to the Monte Carlo approximation, which is undesirable. On the other hand, it equalizes the particle weights, which increases the accuracy of subsequent importance sampling approximations. The decision of whether or not to resample is typically based on a threshold rule for the variance of the particle weights which can be transformed into an effective particle sample size (ESS):
If the particles have equal weights, then $\widehat{ESS}_n = N$. If one particle has weight $N$ and all other particles have weight 0, then $\widehat{ESS}_n = 1$. These are the upper and lower bounds for the effective sample size. To balance the trade-off between adding noise and equalizing particle weights, we execute the resampling step if $\widehat{ESS}_n$ falls below $N/2$ using a systematic resampling algorithm.
The mutation step changes the particle values. In the absence of the mutation step, the particle values would be restricted to the set of values drawn in the initial stage from the prior distribution. This would clearly be inefficient, because the prior distribution is typically a poor proposal distribution for the posterior in an importance sampling algorithm. As the algorithm cycles through the $N_\phi$ stages, the particle values successively adapt to the shape of the posterior distribution. This is the key difference between SMC and classic importance sampling. The transition kernel $K_n(\theta|\tilde{\theta}; \zeta_n)$ is designed to have the following invariance property:
Thus, if $\hat{\theta}_n^i$ is a draw from $\pi_n$, then so is $\theta_n^i$. The mutation step can be implemented by using one or more steps of a MH algorithm. The probability of mutating the particles can be increased by blocking the elements of the parameter vector $\theta$ or by iterating the MH algorithm over multiple steps. The vector $\zeta_n$ summarizes the tuning parameters of the MH algorithm.
Up to now we imposed minimal conditions on the sequence of bridge posterior distributions. To initialize the algorithm, we implicitly required required that it is possible to sample from the initial distribution $\pi_0(\theta)$, which is typically the prior $p(\theta)$ under likelihood or data tempering, and we required that the stage $N_\phi$ posterior is equal to the target posterior distribution: $\pi_{N_\phi}(\theta) = \pi(\theta)$. While previous applications of Algorithm (ref) in econometrics focused on either data or likelihood tempering, the contribution of our paper is to assess the performance of the model tempering approach. Under model tempering the bridge distributions are constructed as follows. Let $M_1$ be the target model with likelihood function $p(Y|\theta,M_1)$ and let $M_0$ be an approximating model with likelihood function $p(Y|\theta,M_0)$. We define the bridge likelihood functions that are used in Steps 2(a) and 2(c) of Algorithm (ref) as:
It can be easily seen that $\pi_{N_\phi}(\theta) = p(\theta|Y,M_1)$, as required, and that the algorithm is initialized with draws from the $M_0$ posterior $\pi_0(\theta) = p(\theta|Y,M_0)$. The intermediate distributions are obtained by shifting the weight gradually from the $M_0$ posterior to the posterior of the target model $M_1$.
Model tempering distinguishes itself from the two most widely-used tempering schemes, likelihood tempering and data tempering, neither of which involve an approximating model $M_0$. Under likelihood tempering (e.g., herbst2014sequential) the stage $n$ posterior is constructed from a tempered version of the full-sample likelihood function: \[ p_n(Y|\theta) = p(Y|\theta,M_1)^{\phi_n}. \] Under data tempering (e.g., DurhamGeweke2011) the bridge distributions are obtained from a fraction of the sample observations $\pi_n(\theta) \propto p(Y_{1:\lfloor \phi_n T \rfloor}|\theta,M_1)p(\theta)$ or, as in Cai {\em et al.} (2021), by gradually shifting the weight from a short-sample likelihood to a full-sample likelihood: \[ p_n(Y|\theta) = p(Y_{1:T}|\theta,M_1)^{\phi_n}p(Y_{1:T_0}|\theta,M_1)^{1-\phi_n}, \quad T_0 < T. \]
Model tempering is a computationally efficient alternative under two conditions. First, the likelihood evaluation of the target model $M_1$ is computationally costly, whereas the likelihood evaluation of the approximating model $M_0$ is, in relative terms, fast. Second, the likelihood functions of the target and the approximating model have to be sufficiently close such that only a modest number of intermediate stages are required to convert draws from the $M_0$ posterior into draws from the $M_1$ posterior. We provide a more detailed discussion in Section (ref) below.
{\bf Adaptive Tempering Schedule.} Under the adaptive tempering schedule used in Cai {\em et al.} (2021) $\phi_n$ is chosen to target a desired level of the ESS defined in ((ref)). Emphasizing the dependence of the incremental weights on the current tempering coefficient $\phi$, write $\tilde{w}_n^i$ in ((ref)) as \[ \tilde{w}^i(\phi) = \frac{p(Y|\theta^i_{n-1},M_1)^\phi p(Y|\theta^i_{n-1},M_0)^{1-\phi}}{p(Y|\theta^i_{n-1},M_1)^{\phi_{n-1}} p(Y|\theta^i_{n-1},M_0)^{1-\phi_{n-1}}} \] and define \[ f(\phi) = \widehat{ESS}_n(\phi) - \alpha \widehat{ESS}^*_{n-1}, \quad 0 < \underline{\alpha} \le \alpha < 1, \] where $\widehat{ESS}^*_{n-1} = \widehat{ESS}_{n-1}$ if the stage $n-1$ selection step (resampling) was executed and $\widehat{ESS}^*_{n-1} = N$ otherwise. Let $\phi_n^*$ satisfy $f(\phi_n^*)=0$ and define $\phi_n = \min \{\phi_n^*,1\}$.
The parameter $\alpha$, to be specified by the user, is the targeted reduction in ESS. It can be shown that for $\phi > \phi_{n-1}$ the ESS satisfies the inequality $\widehat{ESS}_n(\phi) < \widehat{ESS}^*_{n-1}$. Moreover, $\widehat{ESS}_n(\phi)$ is a strictly decreasing function of $\phi$ such that $f(\phi)=0$ has a unique solution. The adaptive algorithm chooses the tempering schedule to control the deterioration of the ESS statistic. The smaller the slope of the function $\widehat{ESS}_n(\phi)$, the larger the increments in the tempering schedule. The number of stages $N_\phi$ is then endogenously determined and is equal to the stage $n$ at which $\phi_n = 1$.
{\bf Model-Specific Parameters.} It might be the case that not all of the parameters that appear in $M_1$ also affect $M_0$, or vice versa. For instance, in one of our illustrations, $M_1$ is a VAR with SV, whereas $M_0$ is a homoskedastic VAR. Thus, the $M_1$ parameter vector contains additional parameters that govern the dynamics of the SV processes. Partition $\theta' = [\theta_c',\theta_0',\theta_1']$, where $\theta_c$ is the vector of common parameters and $\theta_j$ are parameters specific to model $M_j$. The likelihood functions are given by \[ p(Y|\theta,M_j) = p(Y|\theta_c,\theta_j,M_j), \quad j=0,1 \] and
Consider stage $n=0$ with $\phi_0=0$. Because $\theta_1$ does not enter the $M_0$ likelihood function, its distribution does not get updated in view of the data $Y$ and we can factorize the $M_0$ posterior as follows. \[ \pi_0(\theta) = p(\theta|Y,M_0) = p(\theta_c,\theta_0|Y,M_0) p(\theta_1). \] Thus, the model tempering SMC algorithm starts from posterior draws of $(\theta_c,\theta_0)$ and prior draws from $\theta_1$. The use of prior draws for $\theta_1$ in the absence from any information through $M_0$ is both natural and desirable.
At stage $n=N_\phi$ the SMC algorithm approximates the $M_1$ posterior which, on the enlarged parameter space, is given by \[ \pi_{N_\phi}(\theta) = p(\theta|Y,M_1) = p(\theta_c,\theta_1|Y,M_1) p(\theta_0). \] The ultimate object of interest is, in slight abuse of notation, the marginal posterior \[ p(\theta_c,\theta_1|Y,M_1)=\int \pi_{N_\phi}(\theta_c,\theta_0,\theta_1) d\theta_0 . \] While the SMC sampler generates draws from the joint posterior of $(\theta_c,\theta_0,\theta_1)$, draws from the marginal posterior can be obtained by simply dropping the $\theta_0$ draws. A potential disadvantage of including $\theta_0$ into the definition of $\theta$ is that the SMC algorithm has to turn $\theta_0$ draws from a potentially highly concentrated posterior $p(\theta_0|Y,M_0)$ into draws from a more diffuse prior $p(\theta_0)$, which may require an undesirably large number of steps. Thus, we recommend to simply fix $\theta_0$ at a reasonable value, e.g., the posterior mean or mode from a preliminary estimation of $M_0$, and then drop it from the definition of $\theta$.
{\bf Marginal Data Density Ratio.} The SMC algorithm produces as a by-product an approximation of the marginal likelihood ratio $p(Y|M_1)/p(Y|M_0)$. Note that
where
In turn, it can be shown that
as the number of particles $N \longrightarrow \infty$; see, for instance, herbst2014sequential.
{\bf Tempered $M_0$ Posterior.} Rather than using the full-information posterior under $M_0$ as the proposal density, one can choose to incorporate only a fraction of the information embodied in the posterior under model $M_0$. Suppose that the draws from $M_0$ are generated through an SMC algorithm with likelihood tempering, which is what we are doing in the illustrations in Sections (ref) and (ref). Then we can define
which leads to the initialization
The density $\pi_0(\theta;\psi_*)$ represents the posterior obtained from the tempered $M_0$ likelihood function. Thus, the posterior sampling for the approximating model is terminated at $\phi_{N_\phi}=\psi_* < 1$ instead of $\phi_{N_\phi}=1$. The advantage of this strategy is that for $\psi_*<1$ the density $\pi_0(\theta;\psi_*)$ is more diffuse than the full $M_0$ posterior and may exhibit a greater overlap with the target posterior in applications in which $M_0$ and $M_1$ posteriors differ substantially. Note that for $\psi_*=0$ the model $M_1$ would be estimated by standard likelihood tempering instead of model tempering.
To formalize the discussion of the computational advantage of model tempering we begin by introducing some additional notation. Let $\tilde{N}_0(\psi_*) = N_\phi^0(\psi_*)+1$ be the number of $M_0$ SMC stages to obtain a particle swarm that approximates $\pi_0(\theta;\psi_*)$ in ((ref)), including the initial stage, which draws from the prior $p(\theta)$. For the subsequent $M_1$ model tempering we write the number of stages as $\tilde{N}_1(\psi_*) = N_\phi^1(\psi_*)+1$, again to emphasize the dependence on $\psi_*$. We regard $\psi_*=0$ as $M_1$ likelihood tempering and adopt the convention that $\tilde{N}_0(0)=0$.
In the typical VAR and DSGE model applications for which the model tempering procedure is developed, the runtime of the SMC algorithm is predominantly determined by the time it takes to evaluate the likelihood function of the underlying models. Let $N_*$ be the number of likelihood evaluations per SMC stage. It is given by $N_* = N \cdot N_{MH} \cdot N_{blocks}$, where $N$ is the number of particles, $N_{MH}$ is the number of Metropolis-Hastings (MH) steps during the mutation phase, and $N_{blocks}$ is the number of parameter blocks used in each MH step. Moreover, let $\tau_j$ be the time it takes to evaluate the likelihood function of model $M_j$, $j=0,1$. Then the total runtime is given by
where $\mathbb{I}\{x > a\}$ is the indicator function that is equal to one if $x>a$ and equal to zero otherwise. Under likelihood tempering, i.e., $\psi_*=0$, the likelihood function of $M_1$ has to be evaluated $N_* \tilde{N}_1(0)$ times and there is no need to evaluate the $M_0$ likelihood. Under model tempering with $\psi_* > 0$, the likelihood function of $M_1$ has to be evaluated $\tilde{N}_1(\psi_*)$ times and the likelihood of $M_0$ needs to be evaluated during the $M_0$ SMC run and the $M_1$ SMC run.
As mentioned in Section (ref), we are concerned with the case in which the use of the (tempered) $M_0$ posterior reduces the number stages for the $M_1$ SMC and the $M_1$ likelihood evaluation is substantially faster than the $M_0$ evaluation: \[ \tilde{N}_1(\psi_*) < \tilde{N}_1(0) \; \mbox{for} \; \psi_* > 0 \quad \mbox{and} \quad \tau_0 < \tau_1. \] In the numerical illustrations in Sections (ref) and (ref) we report runtimes of model tempering relative to likelihood tempering:\footnote{We found that this formula approximates the actual runtime reductions well.}
The first ratio on the right-hand side of ((ref)) captures the effect of reducing the number of $M_1$ SMC stages needed to reach the target posterior by starting from the tempered $M_0$ posterior $\pi_0(\theta;\psi_*)$ instead of the prior $p(\theta)$. It does not depend on the relative runtime of the $M_1$ and $M_0$ likelihood evaluations. The second term captures the relative costs of having to evaluate the $M_0$ likelihood function. If $\tilde{N}_1(\psi_*)+\tilde{N}_0(\psi_*) \approx \mbox{const}$ as a function of $\psi_*$, then the second term generates a level shift of ${\cal R} (\psi_*,\tau_0/\tau_1)$. As the likelihood evaluation of $M_0$ becomes costless relative to the $M_1$ likelihood evaluation,
In the limit, the time it takes to estimate $M_0$ becomes irrelevant and the reduction is purely driven by the reduction in the number of SMC stages resulting from using an initial distribution that is closer to the target posterior.
In Figure (ref) we provide a numerical example for the runtime reduction. In the left panel, we plot functions $\tilde{N}_0(\psi_*)$ and $\tilde{N}_1(\psi_*)$ which are obtained from DGP 1 of the VAR-SV illustration in Section (ref). The functions are evaluated at $\psi_* \in \{0.0, 0.2, 0.4, 0.6, 0.8, 1.0\}$. As $\psi_*$ increases, the number of stages used in the $M_0$ SMC, denoted by $\tilde{N}_0(\psi_*)$, rises, whereas the number stages in the $M_1$, $\tilde{N}_1(\psi_*)$, falls. The total number of stages required to reach the target posterior stays approximately constant.
The right panel of Figure (ref) depicts ${\cal R} (\psi_*,\tau_0/\tau_1)$ for various choices of $\tau_0/\tau_1$. In this example the most runtime drastic reduction occurs by moving from likelihood tempering to $\psi_*=0.2$. For $\psi_* \ge 0.6$ the function is essentially flat. In the VAR illustration $\tau_0/\tau_1=1/9$. We reduce the likelihood-evaluation ratio all the way to 0. The figure indicates that the reduction in the ratio creates a modest downward shift of the level because the sum $\tilde{N}_1(\psi_*)+\tilde{N}_0(\psi_*)$ is fairly insensitive to $\psi_*$.
Thus far, we have provided an {\em ex post} evaluation of computational gains that relied on knowing how the number of SMC stages depends on $\psi$ through the functions $\tilde{N}_1(\psi_*)$ and $\tilde{N}_0(\psi_*)$. To conduct an {\em ex ante} assessment, we recommend the researcher first assesses the times $\tau_j$ it takes to evaluate the likelihood function of the two models. Moreover, we recommend for several values of $\psi_*$ to compute the variance (across $i$) of the importance sampling weights
where the $\theta^i$'s are draws from $\pi_0(\theta) \propto p(Y|\theta,M_0)^{\psi_*} p(\theta)$. If there is a $\psi_*>0$ for which this variance is considerably smaller than for the $\psi_*=0$ (prior) weights, then there is potential for a substantial runtime reduction. We further explore the relationship between importance sample and runtime reductions in the context of the VAR and DSGE illustrations in Sections (ref) and (ref).
In the first numerical illustration, we consider an environment in which we can directly control the discrepancy between the approximate posterior and the target posterior. We examine the performance of the model tempering approach as a function of the discrepancy between the posteriors. Starting points are “posterior” densities $p(\theta|Y,M_0)$ (approximate) and $p(\theta|Y,M_1)$ (target). We assume that $\theta$ is scalar and approximate and target density are both Normal. In particular, we hold the target density fixed at $p(\theta|Y,M_1) \sim N(0,1)$ and consider a family of approximating densities $p(\theta|Y,M_1) \sim N(\mu,\sigma^2)$, where $\mu$ ranges from -3 to 0 in 0.5 increments and $\sigma$ ranges from 0.2 to 2 in 0.2 increments.
Because in this example we do not construct the posterior density explicitly from a prior distribution and a likelihood function, we let\footnote{This is a slight abuse of notation because $p_n(\theta|Y)$ is not a properly normalized density of $\theta$.} \[ p_n(\theta|Y) = p(\theta|Y,M_1)^{\phi_n} p(\theta|Y,M_0)^{1-\phi_n} \] and define the incremental weight $\tilde{w}_n^i$ in ((ref)) as \[ \tilde{w}_n^i = \frac{p_n(\theta_{n-1}^i|Y)}{p_{n-1}(\theta_{n-1}^i|Y)}. \] We run the SMC Algorithm (ref) with $N=1,000$ particles, $\alpha = 0.95$, and use $N_{MH}=1$ iteration of a single-block random walk Metropolis-Hastings (RWMH) algorithm in the mutation step with $c_0=0.5$, targeting an acceptance probability of 0.25. The implementation of the mutation step is described in more detail in the Online Appendix.
Figure (ref) illustrates how the particle swarm moves from an approximate posterior to the target posterior, despite very little overlap between the two densities. We overlay the target posterior density, $N(0,1)$, an approximate density that is used in this example to initialize the algorithm, $N(-3,0.2)$, and a histogram constructed from the stage $n$ particle swarm. For $n=0$ (left panel) the particle swarm represents the approximate posterior $p(\theta|Y,M_0)$, and for $n=N_\phi=71$ (right panel) it represents the target posterior density $p(\theta|Y,M_1)$. In the center panel of the figure we consider the value of $n=40$ for which the particle swarm represents a weighted geometric mean of the two densities with $\phi_{40}=0.9$.
In Figure (ref) we illustrate the performance of the SMC algorithm across $N_{run}=100$ runs for the different choices of the approximate posterior. To graphically present the results, we mapped $(\mu,\sigma)$ into a discrepancy measure. In principle one could use the variance of the importance weights defined in ((ref)). However, it turns out that we consider fairly large discrepancies between $M_0$ and $M_1$ for which the population variance of the importance weights is infinite. Thus, in this section we use an alternative discrepancy measure defined as one minus the area under the minimum of the two densities: \[ {\cal D}(M_0,M_1) = 1- \int \min \, \big\{ p(\theta|Y,M_0),\, p(\theta|Y,M_1) \big\} d\theta. \] By construction $0 \le {\cal D}(M_0,M_1) \le 1$.
Panel (1,1) shows the average number of stages, $N_\phi$, as a function of ${\cal D}$. In general, the smaller the discrepancy ${\cal D}(M_0,M_1)$, the lower $N_\phi$. Because multiple combinations of $(\mu,\sigma)$ can lead to the same ${\cal D}(M_0,M_1)$, the graph associates multiple $N_\phi$ values with particular values of the discrepancy. For the same level of overlap, approximate densities with a larger standard deviation require fewer stages. This observation provides a justification to start the SMC algorithm from a tempered posterior of the approximating model rather than the full posterior.\footnote{This observation is related to the well-known importance sampling result that the proposal density should have fatter tails than the target density, .e.g., see Geweke1989.}
In Panel (1,2) we depict the tempering schedules for $\mu = -1$ and various values of $\sigma$. Because the mean of the approximating density is different from the mean of the target density, increasing the standard deviation $\sigma$ from 0.2 to 2.0 increases the overlap of the two densities, decreases ${\cal D}(M_0,M_1)$, and leads to a steeper tempering schedule. The runtime pattern in Panel (2,1) mirrors the pattern of the average number of stages, because the runtime increases linearly in the number of stages. Finally, we show the standard deviation of the Monte Carlo approximation of $\mathbb{E}[\theta|Y,M_1]$ as a function of ${\cal D}(M_0,M_1)$ in Panel (2,2). Here no clear relationship with the discrepancy between approximate and target posterior emerges. Because the tempering schedule is chosen adaptively, the accuracy can be as good for large mismatches as it can be for small discrepancies, but it takes more time in the former case.
We deduce from this example that (i) the speed of the model tempering approach depends on the discrepancy between the approximating and the target density. (ii) Except for the additional runtime, the algorithm still works for fairly large discrepancies between the two densities. (iii) It might be desirable to start from a tempered rather than the full posterior of the approximating model.
This section demonstrates the benefits of model tempering in the context of a VAR with SV. The illustration is based on the VAR analysis in AruobaMlikotaSchorfheideVillalvazo2021, except that we do not include a censored endogenous variable. The VAR model that is used as data generating process (DGP) and then estimated based on simulated data is presented in Section (ref). The parameterization of the VAR and the tuning of the SMC algorithm are summarized in Section (ref). The numerical results are discussed in Section (ref).
Model $M_1$ is taken to be a bivariate VAR(1) with stochastic volatility:
where $\Sigma$ is a symmetric positive definite matrix and $\mbox{chol}(\cdot)$ is the lower-triangular Cholesky factor. Let $d_t = [d_{1,t},d_{2,t}]'$ and assume that its elements evolve according to
The presence of stochastic volatility renders this model nonlinear. However, the conditional linearity makes the likelihood evaluation relatively straightforward. We use a Bootstrap Particle Filter (BSPF) with $M_{bspf}=100$ particles, as outlined in the Online Appendix, to sequentially integrate out the latent volatility states. The BSPF likelihood evaluation can be conveniently integrated into the SMC sampler described in Algorithm (ref).\footnote{The use of a particle filter to evaluate the likelihood in the SMC posterior sampler results in a SMC$^2$ algorithm, as discussed in ChopinJacobPapas2012.} Moreover, this computational strategy is very similar to a Bayesian estimation approach widely-used for the estimation of nonlinear DSGE models. We will use the BSPF also in Section (ref).
The approximating model $M_0$ is identical to $M_1$ except that it ignores stochastic volatility. It is given by the homoskedastic VAR
In other words, $M_0$ is obtained by setting $\xi_i = 0 \; \forall \; i$. This restriction renders $(\rho_i,\xi_i)$ non-identified. One obtains the standard analytical expression for the likelihood of a VAR, which as a result can be evaluated instantaneously. Thus, one important condition that makes model tempering attractive is satisfied: the evaluation of the likelihood function for the approximating model is considerably faster than the evaluation of the target model's likelihood.
We use a version of the Minnesota prior for $(\Phi_1,\Phi_c,\Sigma)$. The (marginal) prior for each $\rho_i$ is a Uniform distribution, while that for $\xi_i$ is an inverse Gamma distribution. Further details on the prior are provided in the Online Appendix. The prior specification is the same for all DGP parameterizations.
Estimation is conducted on $T=100$ observations simulated from model $M_1$. We use the following parameterization for $(\Phi_1,\Phi_c,\Sigma)$: \[ \Phi_1 =
, \quad \Phi_c =
, \quad \Sigma =
\cdot
=
. \] The closeness of the posteriors under the target model (VAR with SV) and the approximating model (homoskedastic VAR) depends on the parameterization of the stochastic volatility processes. We consider three different parameterizations which are summarized in Table (ref). Under DGP 1 (baseline) the standard deviations of the log volatility innovations are relatively small. This implies that the $\ln d_{it}$s only exhibit modest time variation and the homoskedastic specification provides a good approximation. Under DGP 2 the volatility innovations have larger standard deviations but the log volatility processes are less persistent, implying large yet short-lived swings in volatility. Finally, DGP 3 combines the baseline values for $\rho_i$ with the large values of $\xi_i$ also considered under DGP 2, implying the largest distance between the approximating model and the target model. This is confirmed in Figure (ref). The panels in the top row show the volatility paths, $d_{1t}$ and $d_{2t}$, and the bottom row illustrates the resulting discrepancy between the $M_0$ and $M_1$ posteriors, using the parameter $\Phi_{1,21}$ as an example.
We use an adaptive tempering schedule with $\alpha=0.95$, as described in Section (ref), to ensure that the number of SMC stages and hence the computational time adjust endogenously to the distance between the proposal and the target density. We initialize the SMC algorithm based on the following set of tempered $M_0$ posteriors: \[ \pi_0(\theta) \propto p(Y|\theta,M_0)^{\psi_*} p(\theta), \quad \psi_* \in \{0.0,0.2,0.4,0.6,0.8,1.0\}, \] where $\psi_*=0$ corresponds to likelihood tempering, i.e., estimation of $M_1$ without using information from model $M_0$. Higher values for $\psi_*$ increasingly tilt the proposal density away from the prior distribution towards the posterior under the proxy model $M_0$. For $\phi_*=1$, the proposal density coincides with the posterior under $M_0$. This is illustrated in Figure (ref) which shows along with the $M_1$ target posterior the sequence of approximating posterior distributions for $\Phi_{1,21}$ under DGP 2 obtained from the $\psi_*$-tempered $M_0$ likelihood function.
The number of particles in the SMC sampler is set to $N=500$. For each DGP and $\pi_0(\theta)$ we run the SMC algorithm $N_{run}=200$ times. We subsequently report averages across the $N_{run}$ runs and assess the variance of the Monte Carlo approximations across runs.
The left panel of Figure (ref) plots average (across multiple SMC runs) Monte Carlo approximates of the log MDD of model $M_1$, $\ln p(Y|M_1)$, under DGP 1 as a function of $\psi_*$, i.e. as a function of the degree of model tempering used in the construction of $\pi_0(\theta)$. The flat line confirms that the Monte Carlo approximations are the same regardless of $\pi_0(\theta)$, as we are numerically approximating the same object regardless of $\psi_*$. The shaded area is a 95% credible band for the log MDD. Figure (ref) in the Online Appendix confirms that the Monte Carlo approximations for posterior mean, variance, 5th and 95th percentiles of the VAR parameters are also invariant to $\psi_*$. Moreover, the result holds not just under DGP 1, but also the other two DGPs (not shown in the figures).
The right panel of Figure (ref) shows the standard deviation of the Monte Carlo approximation of the log MDD as a function of $\psi_*$ for the three different DGPs. For DGP 1 and DGP 3 the standard deviations are weakly decreasing in $\psi$. The biggest drop occurs between $\psi_*=0$ and $\psi_*=0.2$. For DGP 2 the profile is approximately flat, that is, based on the discrepancy between approximating density and target density, the algorithm adjusts the number of stages to keep the accuracy of the Monte Carlo approximation approximately constant.
We now proceed by examining how $\psi_*$ affects the runtime of the SMC algorithm. The main result is presented in Figure (ref). The left panel compares the runtime profiles (normalized by the LT runtime) across the three DGPs. Incorporating information from model $M_0$ in the construction of the proposal density drastically reduces the runtime. The ratio of likelihood evaluation times for $M_0$ and $M_1$ is $\tau_0/\tau_1 = 9.14$. For DGP 1 the runtime monotonically decreases as the proposal is increasingly tilted towards the posterior of the approximate model $M_0$, with the largest reduction of close to 80% obtained for $\psi_* \in \{0.6,0.8,1.0\}$. The runtime reduction is largest for DGP 1, followed by DGP 2, while DGP 3 is associated with the smallest reduction.
The steepest decrease in runtime occurs at $\psi_*=0.2$, which means that incorporating just a bit of information from the approximate model $M_0$ when constructing the proposal density for the posterior of model $M_1$ goes a long way in reducing the runtime. Adding more information helps little at best and might be even counterproductive, as is the case for DGP 3. In principle, for an even larger distance between the posteriors under the two models, it is conceivable that there is no runtime reduction at all. In this case, the second condition stated in Section (ref) would be violated.
The runtime benefits of model tempering decrease with the distance between the posteriors under the target model $M_1$ and the approximate model $M_0$.\footnote{The runtime for DGP 3 also exhibits the most variation across multiple runs; see Figure (ref) in the Online Appendix. Depending on the run, there could be many or only very few particles in the small area to which both posteriors assign some positive probability mass.} In our VAR application, this distance increases with the nonlinearities generated by the SV specification, which are strongest for DGP 3. This distance is visualized in Figures (ref), (ref) and (ref) in the Online Appendix. While for DGP 1 all marginal posteriors align very well for the two models, for DGP 3 there are some parameters for which there is little overlap of probability mass between the two distributions.
We discussed in Section (ref) that this distance could be assessed {\em ex ante}, without having completed the SMC run for $M_1$, by computing the variance of the importance sampling weights defined in ((ref)). The variances as a function of $\psi_*$ for the three DGPs are depicted in the right panel of Figure (ref). The variance profiles look very similar to the runtime profiles in the left panel. For $\psi_*=0$ the variance is approximately equal to the number of particles minus one, $N-1$, which means that one of the particles has weight $N$ and the remaining particles have weight zero. For $\psi_*=0.2$ the variance is considerably lower: it is 75 for DGP 1 and 357 for DGP 3.
To summarize, in this VAR-SV application model tempering is able to reduce the relative runtime by between 27% (DGP 3) and 79% (DGP 1) and increase the precision of Monte Carlo approximations (DGP 1 and 3).
Finally, we consider the estimation of a nonlinear DSGE model, which is computationally demanding for two reasons: first, the model needs to be solved and, second, the evaluation of the likelihood function requires a nonlinear filter. For the latter task, we will use a particle filter similar to the one used to estimate the VAR with SV in Section (ref). In our numerical illustration, we consider a real business cycle model with asymmetric quadratic capital adjustment costs. The adjustment cost parameters let us control the degree of nonlinearity. The model economy is described in Section (ref), the configuration of the simulation experiment is summarized in Section (ref), and the simulation results are presented in Section (ref).
The model economy consists of a representative household and a representative firm. The household consumes $C$, supplies labor in the amount of $L$, and owns the capital stock $K$. The firm hires labor and rents capital to produce a single good that can be used for consumption and investment. The model dynamics can be described as the solution to the following social planner problem:
Households derive utility from consumption and disutility from labor. The parameter $\beta$ is the discount factor, $\tau$ determines the risk aversion, and $\nu$ is the Frisch labor supply elasticity. The parameter $\alpha$ is the capital share parameter, and $\delta$ the depreciation rate. Total factor productivity $Z$ and the preference process $B$ evolve exogenously according AR(1) laws of motion:
Thus, $\varepsilon_z'$ can be thought of as a supply and $\varepsilon_x'$ as a demand shock. Here we adopt the convention that for a variable $X$, $X_*$ denotes the steady state and $\hat{x}$ denote log deviations from the steady state.
For the adjustment cost function we use a linex function which is asymmetric:
The parameter $\phi_1$ controls the overall level of adjustment costs and $\phi_2$ determines the asymmetry. Notice that as $\phi_2 \longrightarrow 0$ the adjustment costs become quadratic around the replacement investment level at which $K'/K=1$. If $\phi_2 > 0$, then it is more costly to reduce the capital stock than it is to augment the capital stock.
Model $M_1$ refers to a nonlinear solution of the RBC growth model, obtained using a second-order perturbation around the steady state. The approximate model $M_0$ is obtained by conducting a first-order linearization. The models are estimated based on observations for output, investment, and hours worked. We denote the observed variables by an o-superscript. The measurement equations, now with $t$ subscripts, take the form:
where $\ln Y_t$, $\ln I_t$, and $\ln L_t$ are the model-implied series and the $\eta_t$s are measurement errors. The measurement errors facilitate the use of a particle filter in combination with the nonlinear DSGE model solution. Moreover, they help to overcome the singularity problem generated by fitting a DSGE model with two shocks to three observables. We include the measurement errors in both data generation and estimation and fix their standard deviations such that the variance of the measurement error is approximately 5% of the variation of the series $\ln Y_t$, $\ln I_t$, and $\ln L_t$, respectively.\footnote{The values that we use are $\sigma_Y=.006$, $\sigma_I = .004$, and $\sigma_L = .004$.}
To facilitate the estimation, we reparameterize the model as follows. First, we express the discount factor $\beta$ as a function of an annualized real interest rate (in percentages) $r = 400(1/\beta-1)$. Second, instead of parameterizing the model in terms of steady states of the exogenous processes $(Z_*,B_*)$, we use the steady states of output and labor, $(Y_*,L_*)$, which we set equal to one for both data generation and estimation. Because our observations $\ln Y_t^o$, $\ln I_t^o$, and $\ln L_t^o$ do not contain direct information on the steady state interest rate and the amount of investment necessary to replace depreciating capital stock, we fix $r$ and $\delta$ at their true values. We collect the parameters that are being estimated in the vector $\theta$: \[ \theta = [\tau, \nu,\alpha,\phi_1,\phi_2,\rho_z,\rho_b,100 \sigma_z,100 \sigma_b]'. \]
As in Section (ref), the estimation is conducted using data simulated from model $M_1$. The parameterization of the DGP is summarized in Table (ref). Most of the parameter values are similar to values commonly found in the DSGE model literature, except that we scale up the shock standard deviations and use fairly large asymmetric adjustment costs by setting $\phi_1 = 50$ and $\phi_2 = 200$. We plot simulated sample paths in the Online Appendix. The length of the estimation sample is $T=80$. The remaining columns of Table (ref) describe the prior distribution for the Bayesian estimation.
In the SMC algorithm we use $N=1,000$ particles to represent the distribution of $\theta$, $N_{MH}=2$ Metropolis-Hastings steps in the mutation, and an adaptive tempering schedule with $\alpha=0.95$. While the likelihood function associated with $M_0$ can be evaluated with the Kalman filter, a nonlinear filter is required to compute the likelihood function of $M_1$. We use the same BSPF that was used in Section (ref) for the VAR estimation with $M_{bspf}=2,000$ particles.
As before, we consider $\psi_* \in \{0.0, 0.2, 0.4, 0.6, 0.8, 1.0\}$, where $\psi_* = 0$ corresponds to $M_1$ likelihood tempering. The relative time it takes to evaluate the $M_0$ and $M_1$ likelihood functions is $\tau_0/\tau_1 \approx 1/ 109$. This ratio depends on the number of particles $M_{bspf}$ used in the BSPF. Doubling $M_{bspf}$ would approximately double $\tau_1$ and cut the ratio in half. The numerical results from a single $N_{run}=1$ run of the model tempering SMC algorithm for the various values of $\psi_*$ are presented in Figure (ref). The top left (1,1) panel shows the log MDD approximation, which is approximately constant as a function of $\psi_*$. This plot confirms that regardless of $\psi_*$ the SMC algorithm delivers the same approximations of the posterior distribution.
The tempering schedules are plotted in Panel (1,2). Starting from a (tempered) $M_0$ posterior drastically reduces the number of stages needed to reach the target posterior. This is consistent with the information provided by the importance weight variance in Panel (2,2). Reweighting draws from the prior $p(\theta)$ ($\psi_*=0$) to target the $M_1$ posterior would lead to a degenerate distribution of weights, whereby the weight of one draw is equal to one and all other weights are equal to zero. Starting with draws from the $\psi_*=0.2$ tempered $M_0$ posterior reduces the importance weight variance from $N-1=999$ to 37. Raising $\psi_*$ toward one, lowers the variance further.
Panel (2,1) depicts the runtime of the model-tempered SMC relative to the $M_1$ likelihood-tempered SMC. The biggest drop, from 1.0 to 0.3 occurs by raising $\psi_*$ from 0.0 to 0.2. Subsequent gains are smaller and the curve essentially turns flat from 0.8 onwards, converging to 0.12. In terms of absolute runtimes, on a Windows workstation with Intel(R) Xeon(R) CPU E5-2687 at 3.10GHz using 8 out of 10 cores model tempering with $\psi_*=1$ reduces the runtime of our JULIA code from 655 to 80 minutes. In the Online Appendix we plot the target and approximate (marginal) posterior densities for the DSGE model parameters. Despite the large $\phi_1$ and $\phi_2$ values, the nonlinearity generated by the DSGE model is not particularly strong and $M_0$ and $M_1$ posterior distributions are quite similar. Thus, when starting from $\psi_*=0.8$ or $\psi_*=1$ only minimal adjustments are required to turn the $M_0$ posterior draws into $M_1$ posterior draws which leads to substantial computational gains.
The implementation of posterior samplers for Bayesian inference often requires the explicit evaluation of likelihood functions. Likelihood calculations for macroeconometric models can be computationally demanding, because it may take a long time to solve the underlying structural model or it may be time-consuming to integrate out latent state variables. In this paper we documented how an SMC algorithm with model tempering can speed up posterior sampling for a VAR with stochastic volatility and a nonlinear DSGE model. The method is suitable for applications in which the likelihood evaluation for the target model is computationally costly and there is an approximating model for which the likelihood evaluation is fast and that generates a posterior that is not too different from the posterior of the target model.
\setcounter{page}{1}