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.
120,611 characters · 18 sections · 54 citation commands
Efficient Likelihood-based Estimation via Annealing for Dynamic Structural Macrofinance Models
In macroeconomics and finance, dynamic structural models provide a convenient theoretical framework to explain economic fluctuations, government policies, and asset prices movements.\footnote{For example, dynamic stochastic general equilibrium (DSGE) models, spanning from real business cycle models rbc1982 to new Keynesian models \citep*[see, e.g.,][]{christiano2005, smets2003}, are commonly used to explain and predict comovements of aggregate macroeconomic fundamentals over the business cycle; consumption-based asset pricing models seek to relate agents' consumption to asset prices and explore fundamental determinants of asset prices lucas1978.} A common feature of these models is that agents' decision rules are derived from assumptions of preferences and technologies by solving intertemporal optimization problems. In empirical applications, these models need to be numerically solved and estimated on real macroeconomic and financial data. For a long time, the literature has been following an approximate approach that these models are first log-linearized and cast into approximate linear state-space models and are then estimated using either likelihood-based methods or moment-based methods.\footnote{The estimation of structural models in macroeconomics and finance follow very different directions. In macroeconomics, the log-linearized models are usually estimated via Bayesian MCMC methods; see, for example, \citet*{herbst2016} and \citet*{schorfheide2016chapter}. In finance, the typical practice involves formulating a set of first-order optimality conditions that map onto moment-based estimation of the key parameters of interest; see, for example, \citet*{bansal2007res}, \citet*{bansal2016jme}, and \citet*{gallant2019rfs}.}
It is now known that second-order approximation errors in linear solutions of dynamic economic models have first-order effects on likelihood estimation \citep*{fernandez2006ecta, fernandez2007estimating}, and that errors in the likelihood estimation are compounded with the sample size. In recent work, \citet*{pohl2018higher} show that log-linearization of consumption-based asset pricing models with long-run risks may result in economically significant errors when state variables are persistent. However, when models are solved using more accurate numerical methods such as high-order perturbation methods \citep*{schmitt2004} or projection methods judd1992projection, the solved models become non-linear and/or non-Gaussian state-space models with high-dimensional and complex structures. The aim of this article is to propose an efficient econometric toolbox based on sequential Monte Carlo (SMC) methods that facilitate likelihood-based estimation when non-linear numerical methods are used to solve models.
Our contribution is threefold. First, we propose a novel methodology to efficiently estimate the likelihood function by building on state-of-the-art SMC methods, also known as particle filters. While particle filters have been applied to estimate the likelihood of dynamic economic models fernandez2007estimating, deJong2013Restud, herbst2019tempered, fulop2020bayesian, the successful application of this approach in model estimation remains challenging due to the large variance of the resulting likelihood estimator. A recent work by heng2020controlled aims to address this issue using ideas from the optimal control literature. Their proposed methodology, termed as controlled SMC, constructs globally optimal proposal distributions that take the {entire sequence of observations} into account using approximate dynamic programming schemes. However, for dynamic structural macrofinance models with complex structures and highly informative observations, controlled SMC can easily fail due to serious numerical instabilities, particularly when the initial proposal distributions are far from optimality.
We develop a new methodology, termed as annealed controlled SMC, which prevents such numerical issues from manifesting. The central idea is to rely on annealing to gradually introduce information from the observations and iteratively refine the proposal distributions as the inverse temperature increases. The original controlled SMC method of heng2020controlled can be seen as a special case of our approach with the inverse temperature fixed at one. Our proposed methodology yields a zero-variance likelihood estimator when the optimal proposal distributions are attainable. Practical implementation requires function approximations, hence resulting in sub-optimal proposal distributions. We provide a theoretical analysis to characterize the quality of our proposal distributions, which elucidates the properties and advantages of annealed controlled SMC. Our results reveal the importance of adopting the annealing procedure. The use of annealing has been explored in online settings to construct locally optimal proposal distributions that take the next observation into account herbst2019tempered,godsill2001improvement. In an offline context where globally optimal proposals are desired, the use of annealing to alleviate numerical instabilities was mentioned in earlier work by scharth2016particle, but the main and crucial difference is that their proposal learning procedure does not exploit the proposals learnt at lower inverse temperatures.
Second, to perform full Bayesian inference, we employ annealed controlled SMC within a novel adaptive SMC$^2$ algorithm that recursively approximates a sequence of annealed posterior distributions of the latent states and model parameters. Compared to existing annealed SMC$^2$ routines such as \citet*{duan2015density} and \citet*{svensson2018learning}, our adaptive approach has two key features. The flexibility of our framework allows us to utilize computationally inexpensive SMC methods (e.g., bootstrap particle filter) at low inverse temperatures, thereby efficiently ruling out unlikely parameters at early stages of the algorithm, and switch to annealed controlled SMC method, which is more computationally demanding but accurate, at higher inverse temperatures when the algorithm considers more promising regions of the parameter space. Moreover, under suitable assumptions, we establish law of large numbers and central limit theorems for our estimators of posterior expectations and the model evidence as the number of parameter particles goes to infinity, for any given number of state particles larger than one.
Compared to MCMC methods or particle MCMC methods \citep*{andrieu2010particle} that sample from the posterior distribution using a single Markov chain, our algorithm shares the many advantages of SMC samplers del2006sequential,dai2020invitation. This includes parallelism over the parameter particles, automated tuning of the inverse temperatures and proposal transitions in the parameter space, and an estimator of the model evidence that facilitiates model comparison.
Third, we apply our adaptive SMC$^2$ with annealed controlled SMC nested to estimate two popular dynamic structural models in macroeconomics and finance. The first model is a prototypical new Keynesian dynamic stochastic general equilibrium model (DSGE) that has been studied in woodford2003, an2007bayesian and herbst2016; the second is a consumption-based long-run risk asset pricing model, in which consumption volatility is modelled as an autoregressive gamma process gourieroux2006autoregressive instead of an autoregressive process that is typically adopted in the standard long-run risk model \citep*[see, e.g.,][]{bansal2004risks, bansal2012empirical}. The first application is described in the main article and the second in the Appendix for the sake of brevity.
We conduct extensive simulation studies based on the above new Keynesian DSGE model that is solved using the log-linearized method sims2002solving or the second-order perturbation method with pruning \citep*{schmitt2004, schmitt2007optimal, fernandez2018res}. For the linearized model, whose likelihood function can be evaluated exactly using a Kalman filter, we find that annealed controlled SMC with moderate number of particles delivers a likelihood estimator with negligible variance. In contrast to the bootstrap particle filter gordon1993novel, the performance of annealed controlled SMC is very robust to the magnitude of the standard deviation of the measurement error. For the non-linear model, whose likelihood function is intractable, we find that the variance of annealed controlled SMC log-likelihood estimator with modest number of particles is about four orders of magnitude smaller than that of the bootstrap particle filter with very large number of particles, and that the magnitude of variance reduction holds even for very small standard deviation of the measurement error. Furthermore, we show that the likelihood estimate from the approximate log-linearized model is much smaller than that of the desired non-linear model obtained with annealed controlled SMC.
We then employ our proposed methodology to estimate the new Keynesian DSGE model on real macroeconomic and financial data. While the log-linearized approximation of this model has been investigated in an2007bayesian and empirically estimated in herbst2016, the corresponding non-linear model has not been fully studied yet in the literature.\footnote{In a recent paper by aruoba2020wp, a variant of this new Keynesian DSGE model is solved using a piecewise-linear approximation method and is estimated using a particle MCMC method.} We estimate the resulting non-linear model under the second-order perturbation method with pruning based on three observables: quarterly per capita GDP growth rate, quarterly inflation, and quarterly annualized interest rate. To avoid the issue of having a zero lower bound on the interest rate, we focus on the pre-crisis sample, ranging from 1983:Q1 to 2007:Q4 with a total of 100 observations.
We estimate the model in two settings: when the standard deviations of the measurement errors are fixed at 20% of the corresponding sample standard deviations, as is commonly considered in the literature, and when they are treated as free parameters to be inferred using the data. We notice that while the estimates suggest a low degree of price rigidity when the standard deviations of the measurement errors are fixed, the opposite conclusion is reached when they are treated as free parameters. This suggests the practice of fixing the standard deviations of the measurement errors may distort model implications. Furthermore, we find that the estimated standard deviations are very different from the fixed values and suggest that the model fits the interest rate better than the output growth rate and the inflation rate.
The remainder of the paper is organized as follows. Section (ref) describes the models of interest and its state-space representation. Section (ref) introduces our annealed controlled SMC method to efficiently estimate the likelihood of these models. Section (ref) details our adaptive SMC$^2$ algorithm for parameter inference. Section (ref) presents an application on a non-linear new Keynesian DSGE model; another application on a non-linear and non-Gaussian consumption-based long-run risk asset pricing model is given in the Appendix. Finally, Section (ref) concludes the paper. Supplementary details and proofs of all theoretical results are provided in the Appendix.
In most dynamic structural macrofinance models, the agents' decision rules are derived from assumptions of preferences and technologies by solving intertemporal optimization problems. For estimation and empirical applications, these models are first numerically solved and then cast into the framework of state-space models. For time $t=1,\ldots,T$, let $s_t\in\mathbb{S}\subseteq\mathbb{R}^d$ denote a vector of latent state variables, which may include both endogenous and exogenous variables, and $y_t\in\mathbb{Y}\subseteq\mathbb{R}^{d_y}$ denote a vector of observations. Given a vector of unknown parameters $\theta\in\Theta\subseteq\mathbb{R}^{d_{\theta}}$, solving a dynamic structural model gives the following relations
where $\Phi_{\theta}^{(0)}$, $\Phi_{\theta}$ and $\Psi_{\theta}$ are non-linear functions from the model solution, $\varepsilon_t$ and $u_t$ represent exogenous shocks and observation noise, respectively. When model-implied observables are deterministically related to state variables, $u_t$ usually captures the measurement error of observations.
We view the sequence of states $(s_t)_{t=0}^T$ as a Markov chain on the state-space $\mathbb{S}$ evolving according to
where $\mu_{\theta}(ds_0)$ denotes an initial distribution and $f_{\theta}(ds_t | s_{t-1})$ a Markov transition kernel on $\mathbb{S}$ that depend on parameters $\theta$. For most models in the literature, these measures either admit densities with respect to the Lebesgue measure on $\mathbb{R}^d$, or are supported on a lower-dimensional subspace when there are both endogenous and exogenous state variables.
The sequence of observations $(y_t)_{t=0}^T$ are assumed to be conditionally independent given the latent process $(s_t)_{t=0}^T$ and are distributed according to
where $g_{\theta}(y_t |s_{t-1},s_t)$ denotes the observation density that is also parameter-dependent. The dependence on the latent states at the previous and current times in Equation (ref) is convenient when modeling asset returns data; longer time dependencies can also be accommodated in our proposed methodology.
Given an observation sequence $y_{1:T}=(y_t)_{t=1}^T\in\mathbb{Y}^{T}$, the complete likelihood is {
}To perform parameter inference, we have to integrate out the latent process to compute the likelihood function
and to estimate the latent states, we need to characterize the smoothing distribution
To facilitate computation and prevent numerical instabilities that often arise when working with dynamic structural macrofinance models, it will be beneficial to gradually introduce the influence of the observations $y_{1:T}$. We do so by defining, for an inverse temperature $\lambda\in[0,1]$, the annealed distribution
and the corresponding likelihood and smoothing distribution as {
}
By gradually increasing $\lambda$, we introduce a path of distributions between the law of the latent process $p(ds_{0:T}|\theta)$ and the desired smoothing distribution $p(ds_{0:T} | y_{1:T}, \theta)$ in Equation (ref). In the context of state-space models, similar motivations can be found in \citet*{godsill2001improvement}, \citet*{svensson2018learning}, and \citet*{herbst2019tempered}. We note that the likelihood $p(y_{1:T} | \theta, \lambda)$ in Equation (ref) is different but related to the tempered likelihood $p(y_{1:T} | \theta)^{\lambda}$, considered in \citet*{duan2015density}, via Jensen's inequality.
For most models of practical interest, the quantities in Equation (ref) are intractable and we have to rely on Monte Carlo approximations. Sequential Monte Carlo (SMC) methods, also known as particle filters, can provide state-of-the-art approximations by simulating an interacting particle system of size $N\in\mathbb{N}$ \citep*{doucet2001introduction,chopin2020introduction}. In what follows, we present a generic description of SMC that encompasses several algorithms in a common framework.
At the initial time, we sample $N$ independent states $(s_0^{(n)})_{n=1}^N$ from a proposal distribution $q_0(ds_0|\theta,\lambda)$ on $\mathbb{S}$. The states are then assigned normalized weights $(W_0^{(n)})_{n=1}^N$ that sum to one according to a weight function $W_0^{(n)}\propto w_0(s_0^{(n)};\theta,\lambda)$. We then sample from the weighted particle approximation $\sum_{n=1}^NW_0^{(n)}\delta_{s_0^{(n)}}(ds_0)$ to multiply states with high weights and discard states that are unlikely. This operation, known as resampling, can be seen as sampling ancestor indexes $(a_0^{(n)} )_{n=1}^N$ from a distribution $r(\cdot|W_0^{(1)},\ldots,W_0^{(N)})$ on $\{1,\ldots,N\}^N$. We will consider multinomial resampling where $(a_0^{(n)} )_{n=1}^N$ are independent samples from the categorical distribution on $\{1,\ldots,N\}$ with probabilities $(W_0^{(n)})_{n=1}^N$. For time $t=1,\ldots,T$, we move the resampled states using a proposal transition kernel $q_t(ds_t|s_{t-1},\theta,\lambda)$ on $\mathbb{S}$, i.e. sample $s_t^{(n)}\sim q_t(\cdot|s_{t-1}^{a_{t-1}^{(n)}},\theta,\lambda)$ independently for $n=1,\ldots,N$. These samples are then assigned normalized weights $(W_t^{(n)})_{n=1}^N$ according to the weight function $W_t^{(n)}\propto w_t(s_{t-1}^{a_{t-1}^{(n)}},s_t^{(n)};\theta,\lambda)$. If $t<T$, we perform resampling by drawing the ancestor indexes $(a_t^{(n)})_{n=1}^N$ from $r(\cdot|W_t^{(1)},\ldots,W_t^{(N)})$. We assume that the smoothing distribution $p(ds_{0:T}|y_{1:T},\theta,\lambda)$ at inverse temperature $\lambda\in[0,1]$ is absolutely continuous with respect to the law of the proposal process
with density written as $p(s_{0:T}|y_{1:T},\theta,\lambda)/q(s_{0:T}|\theta,\lambda)$, and that the choice of weight functions satisfy
The complexity of SMC methods is $\mathcal{O}(NT)$, and its memory requirement is also $\mathcal{O}(NT)$ if all states $(s_t^{(n)})_{t=0,n=1}^{T,N}$ and ancestor indexes $(a_t^{(n)})_{t=0,n=1}^{T-1,N}$ are stored. The latter can be lowered to $\mathcal{O}(T+N\log N)$ using efficient implementations \citep*{jacob2015path}. Given the output of SMC, an unbiased estimator of the likelihood $p(y_{1:T}|\theta,\lambda)$ at the inverse temperature $\lambda\in[0,1]$ is
and a weighted particle approximation of the corresponding smoothing distribution is given by
In the above, each trajectory $s_{0:T}^{(n)}$ is formed by tracing the ancestral lineage of $s_T^{(n)}$, i.e. $s_{0:T}^{(n)}=(s_t^{(l_t^{(n)})})_{t=0}^T$ with particle indexes $(l_t^{(n)})_{t=0}^T$ given by the backward recursion $l_T^{(n)}=n$ and $l_t^{(n)} = a_t^{(l_{t+1}^{(n)})}$ for $t=T-1,\ldots,0$. Using the particle approximation in Equation (ref), smoothing expectations $\int_{\mathbb{S}^{T+1}} \varphi(s_{0:T})p(ds_{0:T}|y_{1:T},\theta,\lambda)$ for any integrable function $\varphi:\mathbb{S}^{T+1}\rightarrow\mathbb{R}$ can be approximated by the weighted average $\sum_{n=1}^NW_T^{(n)}\varphi(s_{0:T}^{(n)})$. Convergence properties of these estimators in the limit of the number of particles $N\rightarrow\infty$ are well-understood del2004feynman,chopin2004central. However, a successful implementation of SMC in practice crucially relies on the choice of proposal distributions in Equation (ref). Poor choices will require prohibitively large number of samples to obtain adequate likelihood and state estimators.
In our framework, the bootstrap particle filter (BPF) of \citet*{gordon1993novel} corresponds to having the law of the latent process in Equation (ref) as proposal, i.e. $q_0(ds_0|\theta)=\mu_{\theta}(ds_0)$, $q_t(ds_t|s_{t-1},\theta)=f_{\theta}(ds_t|s_{t-1})$ for $t=1,\ldots,T$, and the weight functions $w_0(s_0)=1$, $w_t(s_{t-1},s_t;\theta,\lambda)=g_{\theta}(y_t|s_{t-1},s_t)^{\lambda}$ for $t=1,\ldots,T$. It is straightforward to verify that these choices satisfy Equation (ref). Although simple to implement, the efficiency of BPF estimators can be particularly poor in practice if the observations are informative, unless the inverse temperature $\lambda$ is small enough to limit the influence of the observations. The fully adapted auxiliary particle filter (APF) introduced by \citet*{pitt1999filtering} and \citet*{carpenter1999improved} can give better performance by constructing proposals that take the current observation into account. This corresponds to selecting the locally optimal proposals $q_0(ds_0|\theta)=\mu_{\theta}(ds_0)$, $q_t(ds_t|s_{t-1},\theta,\lambda)=p(ds_t|s_{t-1},y_t,\theta,\lambda)$ for $t=1,\ldots,T$, and the weight functions $w_0(s_0)=1$, $w_t(s_{t-1};\theta,\lambda)=p(y_t|s_{t-1},\theta,\lambda)$ for $t=1,\ldots,T$. Exact implementation of APF is not feasible for many non-linear and/or non-Gaussian state-space models, such as the ones to be considered in this paper, as the proposal transitions and weight functions are intractable. Various tractable approximations of APF have been considered doucet2000sequential, finke2020limit.
In this paper, we develop a novel methodology that can significantly outperform both BPF and APF, by constructing globally optimal proposal distributions that take the entire observation sequence $y_{1:T}$ into account. Our approach extends the controlled SMC methodology of heng2020controlled that in turn builds upon the works by \citet*{richard2007efficient, scharth2016particle} and \citet*{guarniero2017iterated}.
In what follows, suppose that we have a given SMC method, defined by a specific choice of proposals $(q_t)_{t=0}^T$ and weight functions $(w_t)_{t=0}^T$ satisfying Equation (ref). For notational ease, we will suppress their dependence on the parameters $\theta\in\Theta$ and the inverse temperature $\lambda\in[0,1]$. The key idea is to express the desired smoothing distribution $p(ds_{0:T} | y_{1:T}, \theta, \lambda) = p(ds_{0} | y_{1:T}, \theta, \lambda) \prod_{t=1}^Tp(ds_{t} | s_{t-1}, y_{t:T}, \theta, \lambda)$ as
for $t=1,\ldots,T$, where the sequence of functions $\psi^*=(\psi_t^*)_{t=0}^T$ are defined by the backward recursion
The notation $q_0(\psi_0^*)=\int_{\mathbb{S}}q_0(ds_0)\psi_0^*(s_0)$ denotes the expectation of $\psi_0^*$ with respect to the distribution $q_0$, and $q_t(\psi_t^*|s_{t-1})=\int_{\mathbb{S}}q_t(ds_t|s_{t-1})\psi_t^*(s_{t-1},s_t)$ denotes the conditional expectation of $\psi_t^*$ under the Markov transition kernel $q_t$. In the case of the BPF, given that $\psi^*$ admits the following probabilistic interpretation
it is sometimes referred to as the backward information filter \citep*{briers2010smoothing}.
By obtaining an approximation $\psi=(\psi_t)_{t=0}^T$ of the backward recursion in Equation (ref) using an approximate dynamic programming method that will be discussed in Section (ref), we can construct a new proposal distribution
by mimicking Equation (ref), i.e. define
Following the terminology in heng2020controlled, we will refer to a sequence of non-negative and bounded functions $\psi$ as a policy and $\psi^*$ as the optimal policy. As the choice of policy is specific to the application of interest, we defer these discussions until Section (ref) and assume $\psi$ is such that the proposals $(q_t^{\psi})_{t=0}^T$ in Equation (ref) can be sampled from, and the expectations $q_0(\psi_0)$, $q_t(\psi_t|s_{t-1})$ for $t=1,\ldots,T$ can be evaluated. To employ these newly constructed proposals within SMC, the appropriate weight functions are given by
which satisfy $w_0^{\psi}(s_0)\prod_{t=1}^Tw_t^{\psi}(s_{t-1}, s_t) = p(s_{0:T},y_{1:T}|\theta,\lambda)/q^{\psi}(s_{0:T}|\theta,\lambda)$. An algorithmic description of the resulting SMC method is detailed in Algorithm (ref). To distinguish between the initial SMC method and the new SMC method induced by a policy, we will refer to the former as uncontrolled SMC and the latter as controlled SMC. From the output of Algorithm (ref), we have an unbiased estimator of the likelihood $p(y_{1:T}|\theta,\lambda)$ (Step 3) and an approximate sample from the smoothing distribution $p(ds_{0:T} | y_{1:T}, \theta, \lambda)$ by sampling from the weighted particle approximation in Equation (ref) (or equivalently selecting an ancestral lineage at Step 4). The efficiency of these approximations will ultimately depend on how well the chosen policy $\psi$ approximates the optimal policy $\psi^*$. A more precise characterization of this relationship will be given in Section (ref). We note that the choice $\psi=\psi^*$ is optimal as Algorithm (ref) would yield a zero variance likelihood estimator, for any number of particles $N$, and an exact trajectory from the smoothing distribution.
Suppose we have a policy $\psi=(\psi_t)_{t=0}^T$ that denotes our current approximation of the optimal policy $\psi^*=(\psi_t^*)_{t=0}^T$ defined in Equation (ref). If we define $\phi^*=(\phi_t^*)_{t=0}^T$ using the backward recursion
it can be shown that $\psi^*=\psi\cdot\phi^*=(\psi_t\cdot\phi_t^*)_{t=0}^T$, where $\psi_t\cdot\phi_t^*$ denotes pointwise multiplication of two functions heng2020controlled. This result shows how policy refinements can be performed and identifies $\phi^*$ as the optimal refinement of the current policy $\psi$. The optimal refinement can be viewed as the solution of an associated Kullback--Leibler optimal control problem and Equation (ref) is the corresponding dynamic programming recursion. As noted by heng2020controlled, drawing this connection allows one to exploit approximate dynamic programming (ADP) methods to approximate the optimal refinement. The following is an ADP scheme to approximate $\phi^*$ by combining function approximation and iterating the backward recursion in Equation (ref).
Let $(s_t^{(n)})_{t=0,n=1}^{T,N}$ and $(a_t^{(n)})_{t=0,n=1}^{T-1,N}$ denote the states and ancestors from running controlled SMC with the current policy $\psi$. At the terminal time $T$, we approximate $\phi_T^*= w_T^{\psi}$ by solving the least squares problem
where $\mathbb{F}_T$ is a function class to be specified. By plugging in the approximation $\phi_T\approx\phi_T^*$ in the iterate $\phi_{T-1}^*=w_{T-1}^{\psi}q_T^{\psi}(\phi_T^*)$, we approximate the function $\varphi_{T-1}=w_{T-1}^{\psi}q_T^{\psi}(\phi_T)$ using the least squares problem
where $\mathbb{F}_{T-1}$ is another function class to be chosen. We then proceed in the same manner until the initial time $0$ to obtain a sequence of functions $\phi=(\phi_t)_{t=0}^T$ approximating the optimal refinement $\phi^*$. The refined policy is then given by the update $\psi\cdot\phi=(\psi_t\cdot\phi_t)_{t=0}^T$.
The above ADP scheme is summarized in Algorithm (ref), where we also give alternative expressions involving the proposals $(q_t)_{t=0}^T$ and weight functions $(w_t)_{t=0}^T$ of the uncontrolled SMC, which we will use for our numerical implementation. The cost of running Algorithm (ref) is $\mathcal{O}(T(NC_{\mathrm{evaluate}} + C_{\mathrm{approx}}))$, where $C_{\mathrm{evaluate}}(d)$ is the cost of evaluating one of the weight functions in Equation (ref), and $C_{\mathrm{approx}}(N,d)$ is the cost of each least squares approximation. In the case of linear least squares, $C_{\mathrm{approx}}$ would be linear in $N$. By selecting parametric function classes $(\mathbb{F}_t)_{t=0}^T$ that depend on coefficients $(\beta_t)_{t=0}^T$, only the estimated coefficients parameterizing the policies have to be stored in practice. We will also choose function classes that are closed under multiplication, so that the refined policy $\psi\cdot\phi$ also lie in the same function classes, with coefficients that can be easily updated. The richness of these function classes will determine the quality of $\psi\cdot\phi$ as an approximation of the optimal policy $\psi^*$.
We can then run a controlled SMC (Algorithm (ref)) with proposals defined by the refined policy $\psi\cdot\phi$ to obtain better likelihood and state estimates. Using this SMC output, one could consider another round of policy refinement with ADP (Algorithm (ref)) to produce more efficient SMC estimates. This iterative procedure that alternates between ADP and SMC, with the inverse temperature fixed at the desired level of $\lambda=1$, is studied by heng2020controlled. However, for complex state-space models with strong non-linearities and highly informative observations, typically of dynamic structural macrofinance models, this approach can easily fail when the initial policy is far from optimality as the states used in the ADP approximation would be in the tails of the smoothing distribution in Equation (ref). To prevent such numerical instabilities from manifesting, we propose a novel iterative procedure in the following section that performs policy refinement as the inverse temperature $\lambda$ gradually increases.
Let $(\lambda_i)_{i=0}^I$ denote an increasing inverse temperature schedule with $\lambda_0=0$ and $\lambda_I\leq 1$ that will be pre-specified. Although the inverse temperature of $\lambda=1$ is desired to estimate the quantities in Equations (ref) and (ref), approximating the intermediate quantities in Equation (ref) for $0\leq\lambda<1$ is also useful when we consider parameter inference in Section (ref).
We begin by running the uncontrolled SMC at inverse temperature $\lambda_0=0$, defined by proposals $(q_t)_{t=0}^T$ and weight functions $(w_t)_{t=0}^T$, and initializing the policy $\psi^{(0)}=(\psi_t^{(0)})_{t=0}^T$ as constant one functions, i.e. $\psi_t^{(0)}=1$ for $t=0,\ldots,T$. In the case of the BPF, this simply generates $N$ trajectories from the latent process of Equation (ref) and initializes the policy at optimality for $\lambda_0=0$. Subsequently, for iteration $i=1,\ldots,I$, we construct new proposals for the next inverse temperature $\lambda_i$ by refining the policy $\psi^{(i-1)}$ from the previous inverse temperature $\lambda_{i-1}$ using ADP. More precisely, we employ Algorithm (ref) using $\psi^{(i-1)}$ as the current policy and the previous SMC output at $\lambda_{i-1}$. Writing $\psi^{(i)}$ as the refined policy, we then run a controlled SMC (Algorithm (ref)) at inverse temperature $\lambda_i$ with new proposals defined by $\psi^{(i)}$. We provide an algorithmic summary of the above methodology in Algorithm (ref), which we refer to as annealed controlled SMC (AC-SMC).
We now explain the rationale behind our proposed methodology and how it differs from existing works. Suppose $\psi^{(i-1)}$ is a good approximation of the optimal policy at $\lambda_{i-1}$, which is the case at initialization. The states generated by the resulting controlled SMC at $\lambda_{i-1}$ (Step 2b) will be in regions of high probability mass under the smoothing distribution $p(ds_{0:T} | y_{1:T}, \theta, \lambda_{i-1})$. If the inverse temperature increase is small, these states should also be in high probability regions under $p(ds_{0:T} | y_{1:T}, \theta, \lambda_{i})$ and therefore are good support points to learn the refined policy $\psi^{(i)}$ using the ADP algorithm (Step 2a). The optimal refinement of $\psi^{(i-1)}$ is given by Equation (ref) and can be rewritten as
for $t = T-1,\ldots,1$. In the preceding equations, we reintroduce parameter and temperature dependence for clarity and assume that the choice of initial proposals $(q_t)_{t=0}^T$ are not temperature dependent.
The weight functions $w_0^{\psi^{(i-1)}}(s_0;\theta)$ and $w_t^{\psi^{(i-1)}}(s_{t-1},s_t;\theta,\lambda_{i-1})$ for $t=1,\ldots,T$ are determined by the quality of the previous ADP approximation and is equal to the residual of the least squares approximation at time $t=0,1,\ldots,T$ in the logarithmic scale. These weight functions should therefore be close to constant functions as we have assumed near optimality of $\psi^{(i-1)}$. The ratio of weights $w_t(s_{t-1},s_t;\theta,\lambda_{i})/w_t(s_{t-1},s_t;\theta,\lambda_{i-1})$, which accounts for the increase in inverse temperature, reduces to $g_{\theta}(y_t|s_{t-1},s_t)^{\lambda_i-\lambda_{i-1}}$ in the case of the BPF. If the inverse temperature increment is small, we can expect these ratios to be well-behaved with small fluctuations. By an inductive argument, the conditional expectation $q_{t}^{\psi^{(i-1)}}(\phi_{t}^*|s_{t-1},\theta)$ should also have the same behaviour. Hence we can expect the optimal refinement of $\psi^{(i-1)}$ to be well-approximated by simple functions. Assuming that the chosen function classes are rich enough and the number of samples $N$ is sufficiently large to obtain good least squares estimation, a good approximation of the optimal policy at $\lambda_{i}$ will then ensure stability of the next iteration.
In contrast to the policy refinement procedure in heng2020controlled, which is based on repeated least squares fitting of residuals with the inverse temperature fixed at $\lambda=1$, our methodology can be seen as an extension that allows one to incorporate changes in temperature. The use of annealing to alleviate numerical instabilities was mentioned in earlier work by scharth2016particle, but the main and crucial difference is that their policy learning procedure does not change across iterations to exploit the policies learnt at lower inverse temperatures.
In practical implementations of AC-SMC (Algorithm (ref)), we can monitor the performance of controlled SMC at each inverse temperature (Step 2b) to evaluate the quality of the ADP approximation (Step 2a). Using the above-mentioned relationship between the weight functions of controlled SMC and the residuals in ADP, we can inspect the variance of the SMC weights $(W_t^{(n)})_{n=1}^N$, using for example the effective sample size criterion \citep*{kong1994sequential}, defined at each time $t=0,1,\ldots,T$ as $1/\sum_{n=1}^N(W_t^{(n)})^2$. Using the output of Algorithm (ref), we obtain an unbiased estimator of the likelihood $p(y_{1:T}|\theta,\lambda_I)$, and an approximate sample from the smoothing distribution $p(ds_{0:T}|y_{1:T},\theta,\lambda_I)$.
We analyze the performance of AC-SMC at each iteration, by supposing that we have a current policy $\psi$, and an approximation $\phi$ of the optimal refinement of $\psi$, obtained using ADP (Algorithm (ref)) at inverse temperature $\lambda\in[0,1]$. To measure the effect policy refinement has on controlled SMC (Algorithm (ref)), we consider the Kullback--Leibler (KL) divergence from the proposal distribution $q^{\psi\cdot\phi}(ds_{0:T})$, defined in Equations (ref) and (ref), to the smoothing distribution $p(ds_{0:T}|y_{1:T},\theta,\lambda)$ in Equation (ref), defined as {
}This choice of KL divergence is motivated by the fact that it characterizes the quality of our proposal distribution in the context of importance sampling chatterjee2018sample, i.e., small KL divergence is both sufficient and necessary for good importance sampling approximations.
We first introduce some notation. We denote by $\mu_t^*$ and $\mu_t^{\psi\cdot\phi}$ the marginal distributions of the smoothing and proposal distributions at time $t=0,\ldots,T$, respectively. Let $(q_t^*)_{t=0}^T$ denote the optimal proposals under the optimal policy $\psi^*$, and define for $t=1,\ldots,T$ the joint smoothing distribution $(\mu_{t-1}^*\times q_t^*)(ds_{t-1},ds_t)=\mu_{t-1}^*(ds_{t-1})q_t^*(ds_t|s_{t-1})$ on $\mathbb{S}\times\mathbb{S}$. Note that the above quantities depend on the parameters $\theta\in\Theta$ and inverse temperature $\lambda\in[0,1]$, even though these dependencies are not made explicit for notational simplicity. Our first result provides a decomposition of the KL divergence, in terms of the logarithmic differences between $\phi^*$ and $\phi$ under the marginal distributions $\xi_0^*=q_0^*$, $\xi_0^{\psi\cdot\phi}=q_0^{\psi\cdot\phi}$, and $\xi_t^*=\mu_{t-1}^*\times q_t^*$, $\xi_t^{\psi\cdot\phi}=\mu_{t-1}^*\times q_t^{\psi\cdot\phi}$ for $t=1,\ldots,T$.
In what follows, all proofs are provided in the Appendix (ref).
Next, we will introduce some assumptions which are needed to derive upper bounds of the terms in Equation (ref). The states $(s_t^{(n)})_{t=0,n=1}^{T,N}$ and ancestors $(a_t^{(n)})_{t=0,n=1}^{T-1,N}$ from controlled SMC with current policy $\psi$ define the following empirical measures {
}for $t=1,\ldots,T$, that are used in the least squares approximations in Equations (ref) and (ref). As the number of particles $N\rightarrow\infty$, these measures $(v_t^{\psi,N})_{t=0}^T$ are consistent approximations of $(v_t^{\psi})_{t=0}^T$, defined recursively as
For each time $t=0,\ldots,T$, we define the $L^2$-norm $\| \varphi \|_{L^2(\nu_t^{\psi})}= \nu_t^{\psi}(\varphi^2)^{1/2}$ and the $L^2(\nu_t^{\psi})$ space\footnote{With equivalent classes defined by functions that agree $\nu_t^{\psi}$-almost everywhere.} as the set of measurable functions $\varphi$ with $\| \varphi \|_{L^2(\nu_t^{\psi})}<\infty$. To study ADP (Algorithm (ref)) in the infinite particle regime, we define $L^2$-projections under the function class $\mathbb{F}_t$ and the distribution $\nu_t^{\psi}$
for $\log \varphi\in L^2(\nu_t^{\psi})$. The following assumption concerns our choice of function classes $(\mathbb{F}_t)_{t=0}^T$ within ADP.
Assumption (ref)$(i)$ ensures the existence of a unique projection in Equation (ref), which can be relaxed by letting $L^2$-projections denote the set of minimizers. The residual errors $(e_t^{\psi})_{t=0}^T$ in Assumption (ref)$(ii)$ describes the flexibility of the chosen function classes for our purpose of learning policies. The next assumption pertains to the relationships between the distributions $\xi_t^{*}$, $\xi_t^{\psi\cdot\phi}$ and $\nu_t^{\psi}$. For distributions $\mu$ and $\nu$ defined on a common measurable space, we will write $\mu\ll\nu$ if $\mu$ is absolutely continuous with respect to $\nu$.
Assumption (ref)$(i)$ requires the current policy $\psi$ to induce a reasonably good approximation of the smoothing marginals. As such a condition is unlikely to be satisfied when $\psi$ is given by constant one functions and the inverse temperature $\lambda=1$, this motivates the use of AC-SMC which allows policy refinement as $\lambda$ gradually increases. Assumption (ref)$(ii)$ is a condition on the quality of the proposals $(q_t^{\psi\cdot\phi})_{t=0}^T$ under the refined policy $\psi\cdot\phi$. We now state our main result, which gives recursive bounds of the logarithmic differences appearing in the KL upper bound in Equation (ref).
Equation (ref) shows how the residual errors $(e_t^{\psi})_{t=0}^T$ in the $N=\infty$ limit propagate backward in time. These errors can be small if the function classes are sufficiently rich, in which case, Theorem (ref) and Proposition (ref) would imply good performance of AC-SMC.
In this section, we consider the problem of parameter and state inference in the Bayesian framework. Let $p(d\theta)=p(\theta)d\theta$ denote our prior distribution on the parameter space $\Theta$. We develop a new methodology to approximate the posterior distribution of the parameters and latent states
and the model evidence $p(y_{1:T})=\int_{\Theta}p(d\theta)p(y_{1:T} | \theta)$. Our approach is to employ an SMC sampler del2006sequential to sequentially approximate distributions and their normalizing constants along the path
for $\lambda\in[0,1]$, in which a nested AC-SMC method (Algorithm (ref)) is exploited to approximate the likelihood $p(y_{1:T}|\theta,\lambda)$ and smoothing distribution $p(ds_{0:T} |y_{1:T}, \theta,\lambda)$. As the inverse temperature $\lambda$ increases from zero to one, Equation (ref) bridges between the prior distribution $p(d\theta)p(ds_{0:T}|\theta)$ and the posterior distribution $p(d\theta, ds_{0:T} | y_{1:T})$ in Equation (ref), and the normalization constant $p(y_{1:T}|\lambda)=\int_{\Theta}p(d\theta)p(y_{1:T} | \theta,\lambda)$ goes from one to the model evidence $p(y_{1:T})$.
Our approach falls in the class of SMC$^2$ methods developed by \citet*{fulop2013efficient}, \citet*{chopin2013smc2}, and duan2015density, but differs in important aspects. In Section (ref), we provide a description of our adaptive SMC$^2$ algorithm and discuss how it differs from and relates to the existing literature. Theoretical justifications of the algorithm are presented in Section (ref).
Our proposed SMC$^2$ methodology, detailed in Algorithm (ref), builds on AC-SMC in Algorithm (ref) and a conditional implementation of controlled SMC (Algorithm (ref) in Appendix (ref)).
In Step 1, we initialize the algorithm at inverse temperature $\lambda_0=0$ by sampling $P$ parameters and state trajectories $(\theta^{(p)},s_{0:T}^{(p)})_{p=1}^P$ from the prior distribution $p(d\theta, ds_{0:T} | y_{1:T},\lambda_0)=p(d\theta)p(ds_{0:T}|\theta)$. Subsequently, for iteration $i\geq 1$, we determine the next inverse temperature $\lambda_i\in(0,1]$ in Step 2, so that the previous bridging distribution $p(d\theta, ds_{0:T} | y_{1:T},\lambda_{i-1})$ provides a good importance sampling approximation of the next bridging distribution $p(d\theta, ds_{0:T} | y_{1:T},\lambda_i)$. In this context, we measure the quality of importance sampling approximations using the effective sample size (ESS) criterion
which is based on the unnormalized weights
and the normalized weights $\Omega_i^{(p)}(\lambda)=\omega_i^{(p)}(\lambda)/\sum_{j=1}^P\omega_i^{(j)}(\lambda)$. The ESS lies between $1$ to $P$: the lower bound is attained when one sample has all the normalized weight, while the upper bound is achieved when all samples have uniform weights.
We adapt the next inverse temperature using the following scheme:
with the convention $\inf\emptyset = 1$ to accommodate the terminal iteration, where $\kappa_{\mathrm{ESS}}\in(0,1)$ is a pre-specified threshold that controls the amount of weight degeneracy. We refer readers to dai2020invitation for discussions on the choice of $\kappa_{\mathrm{ESS}}$ and its impact on the number of iterations required to reach the desired inverse temperature of $\lambda=1$. In contrast to the adaptation scheme proposed in \citet*{svensson2018learning}, which requires storing the entire SMC output for each parameter $\theta^{(p)}$, the unnormalized weight in Equation (ref) is considerably simpler as it only depends on a single trajectory. Moreover, as it can be shown that our ESS criterion is a strictly decreasing and continuous function of $\lambda\in(\lambda_{i-1},1]$, Equation (ref) can be implemented using a simple bisection routine. After determining $\lambda_i$, we compute the resulting importance weights (Step 3a) and perform resampling to focus our computational effort on more likely parameters and trajectories (Step 3b). For notational ease, we have used the same notation $(\theta^{(p)},s_{0:T}^{(p)})_{p=1}^P$ to denote the resulting set of samples after any operation.
In Step 3c, we then learn a policy $\psi^{(p)}$ for each resampled parameter $\theta^{(p)}$ to construct good proposal distributions that approximate the smoothing distribution $p(ds_{0:T} | y_{1:T},\theta^{(p)},\lambda_i)$ at inverse temperature $\lambda_i$. This policy learning step is left intentionally general in the algorithm to accommodate various strategies for different applications. A generic and useful recipe is to set $\psi^{(p)}$ as constant one functions for inverse temperatures $\lambda_i$ that are below a small pre-specified level $\lambda_*\in(0,1)$. The rationale here is that at low inverse temperatures, the performance of uncontrolled SMC would be adequate as the influence of the observations are limited. This allows us to efficiently rule out very unlikely parameters at early iterations of the algorithm. For larger inverse temperatures $\lambda_i>\lambda_*$, it is worthwhile spending the computational overhead to learn the optimal policies in more promising regions of the parameter space.
In Step 3d, given the policy $\psi^{(p)}$ for each parameter $\theta^{(p)}$, we then run a conditional implementation of controlled SMC (Algorithm (ref) in Appendix (ref)) to obtain a likelihood estimator $\hat{p}(y_{1:T}|\theta^{(p)},\lambda_i)$ and a new trajectory $s_{0:T}^{(p)}$. The distinctive feature in the conditional implementation is that the input reference trajectory is conditioned to survive all resampling steps \citep*{andrieu2010particle}, which is necessary for the validity of our approach. The use of conditional SMC within SMC$^2$ was also considered by \citet*{chopin2013smc2}, but for the purpose of increasing the number of state particles $N$ as more observations are assimilated. As we are constructing better SMC proposals by learning optimal policies, the choice of $N$ is less crucial in our setting. Compared to the tempered likelihood approach of duan2015density, our proposed methodology has the flexibility to alter the SMC configuration via conditional SMC and has the benefits of annealing. Note that at this stage, we only have to store the resampled parameters $(\theta^{(p)})_{p=1}^P$, updated trajectories $(s_{0:T}^{(p)})_{p=1}^P$, and likelihood estimators $(\hat{p}(y_{1:T}|\theta^{(p)},\lambda_i))_{p=1}^P$.
Compared to having a single Markov chain to sample from the posterior distribution of Equation (ref), the ability to select tuning parameters of a proposal transition kernel $h_i$, based on existing samples that approximate $p(d\theta|y_{1:T},\lambda_i)$ in Step 4, is an advantage of the SMC$^2$ approach. We refer readers to dai2020invitation and references therein for discussions of various adaptation rules. Our numerical experiments employ Gaussian random walk proposals with covariance matrices that are estimated using the current set of parameter particles $(\theta^{(p)})_{p=1}^P$. Step 4 describes $K\in\mathbb{N}$ particle marginal Metropolis--Hastings (PMMH) moves \citep*{andrieu2010particle} to improve the sample diversity of the resampled parameters. The combination of these steps define a Markov transition kernel that has $p(d\theta, ds_{0:T} | y_{1:T},\lambda_i)$ as its invariant distribution. As the efficiency of PMMH moves crucially depends on the variance of the likelihood estimator doucet2015efficient,sherlock2015efficiency, we learn the optimal policy for each proposed parameter $\theta^{(*)}$ (Step 4b), and run controlled SMC with the resulting policy $\psi^{(*)}$ (Step 4c) to obtain a lower variance likelihood estimator $\hat{p}(y_{1:T}|\theta^{(*)},\lambda_i)$ and a trajectory $s_{0:T}^{(*)}$ with a law that is closer to the smoothing distribution $p(ds_{0:T} | y_{1:T}, \theta^{(*)}, \lambda_i)$. In Step 4d, the proposed parameter $\theta^{(*)}$, trajectory $s_{0:T}^{(*)}$, and likelihood estimator $\hat{p}(y_{1:T}|\theta^{(*)},\lambda_i)$ are then accepted according to the Metropolis--Hastings acceptance probability
When the desired inverse temperature of $\lambda=1$ is reached, we terminate the SMC$^2$ algorithm (Step 5). The number of iterations $I$ required to complete the algorithm is random and depends on the choice of ESS threshold $\kappa_{\mathrm{ESS}}\in(0,1)$. Using the parameters and trajectories $(\theta^{(p)},s_{0:T}^{(p)})_{p=1}^P$ outputted by the algorithm, we can approximate posterior expectations $\pi(\varphi)=\int_{\Theta\times\mathbb{S}^{T+1}}\varphi(\theta,s_{0:T})p(d\theta, ds_{0:T} | y_{1:T})$, for any integrable function $\varphi:\Theta\times\mathbb{S}^{T+1}\rightarrow\mathbb{R}$, with the sample average $\hat{\pi}(\varphi)=P^{-1}\sum_{p=1}^P\varphi(\theta^{(p)},s_{0:T}^{(p)})$. As a by-product of the algorithm, we also have an estimator of the model evidence $\hat{p}(y_{1:T})$, computed using the unnormalized weights in Step 6
Since estimators of the model evidence cannot be easily obtained with a single PMMH chain targeting Equation (ref), this is another strength of our SMC$^2$ approach. In the next section, we will establish consistency properties of our estimators of posterior expectations and the model evidence in the limit of the number of parameters particles $P\rightarrow\infty$, for any choice of the number of state particles $N>1$.
This section concerns asymptotic properties of our estimators of expectations under the posterior distribution in Equation (ref) and the model evidence $p(y_{1:T})$. In the following, the notations $\stackrel{p.}{\rightarrow}$ and $\stackrel{d.}{\rightarrow}$ denote convergence in probability and distribution, respectively. We will establish the weak law of large numbers (WLLN)
as the number of parameter particles $P\rightarrow\infty$, for any number of state particles $N>1$. To quantify the rate of convergence, we will also seek the central limit theorems (CLT)
as $P\rightarrow\infty$ for any fixed $N>1$, where we denote the normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$ as $\mathcal{N}(\mu,\Sigma)$ and its density by $x\mapsto\mathcal{N}(x;\mu,\Sigma)$. Although we will not give explicit expressions of the asymptotic variances $\sigma^2(\varphi)$ and $\sigma^2$ due to the complexity of our SMC$^2$ algorithm, it is clear how they can be derived in our proofs, provided in the Appendix (ref).
We first consider Algorithm (ref) without adaptation, i.e., the inverse temperature schedule $(\lambda_i)_{i=0}^I$ and the proposal transition kernels $(h_{i})_{i=1}^{I}$ are pre-specified and not determined on the fly in Steps 2 and 4. This analysis serves to elucidate various aspects of our SMC$^2$ algorithm without the added complication of adaptation.
To study the adaptive SMC$^2$ algorithm, we first note that the adaptation scheme in Step 2 and Equation (ref) should be seen as a finite sample approximation of a limiting deterministic inverse temperature schedule $(\lambda_i^*)$, defined by the following scheme
initialized at $\lambda_0^*=0$, where the above $\chi^2$-divergence
depends on the importance weight $\omega$ defined in Equation (ref). Next, we formalize the adaptive nature of Step 4 using the framework of beskos2016convergence. At iteration $i$, we consider a parametric family of proposal transition kernels $h_i$, indexed by tuning parameters $\xi\in\mathbb{R}^{s}$. The algorithm determines these tuning parameters by computing the sample average $\xi_i=P^{-1}\sum_{p=1}^PS_i(\theta^{(p)})$ of a summary statistic $S_i:\Theta\rightarrow\mathbb{R}^{s}$, based on a current set of parameter particles $(\theta^{(p)})_{p=1}^P$ approximating $p(d\theta|y_{1:T},\lambda_{i-1})$. Let $m_i$ denote the resulting Markov transition kernel on $\Theta\times\mathbb{S}^{T+1}$ by composing Steps 3c, 3d and 4 of Algorithm (ref), which will be shown to have $p(d\theta,ds_{0:T}|y_{1:T},\lambda_{i})$ as its invariant distribution. We also define the non-negative kernel $Q_i(d\tilde{\theta},d\tilde{s}_{0:T}|\theta,s_{0:T})=\omega(\theta,s_{0:T}|\lambda_{i-1},\lambda_i) m_i(d\tilde{\theta},d\tilde{s}_{0:T}|\theta,s_{0:T})$ on $\Theta\times\mathbb{S}^{T+1}$. To stress the dependence of $Q_i$ on the quantities $\zeta_i=(\lambda_{i-1},\lambda_i,\xi_i)$, we will write $Q_{i,\zeta_i}$. Our assumptions will involve the behaviour of the function $\zeta_i\mapsto Q_{i,\zeta_i}$ and its gradient $\zeta_i\mapsto \nabla_{\zeta_i}Q_{i,\zeta_i}$ at $\zeta_i^*=(\lambda_{i-1}^*,\lambda_i^*,\xi_i^*)$, which corresponds to an idealized algorithm where the limiting inverse temperature schedule $(\lambda_i^*)$ is employed, and tuning parameters are determined by the expectation $\xi_i^*=\int_{\Theta}S_i(\theta)p(d\theta|y_{1:T},\lambda_{i-1})$.
The crux of our arguments to establish Theorems (ref) and (ref) is to first cast the rather involved SMC$^2$ algorithm, without and with adaptation, as a particular SMC sampler and adaptive SMC sampler in the frameworks of del2006sequential and beskos2016convergence that operate on an extended space. This enables us to then invoke convergence results for standard SMC methods del2004feynman,chopin2004central and apply them to our SMC$^2$ estimators. Such an approach that exploits the specific properties of Algorithm (ref) is not applicable to the SMC$^2$ method of duan2015density, which is based on the tempered likelihood.
Dynamic stochastic general equilibrium (DSGE) models have been widely used in macroeconomic research and in central banks for forecasting and policy-making. In this section, we consider a prototypical New Keynesian DSGE model that has been studied in woodford2003, an2007bayesian, and herbst2016. This model has become a benchmark specification for the analysis of monetary policy. Variants of this model have also been studied in \citet*{aruoba2018res} and aruoba2020wp.
In the Appendix (ref), we provide another application to estimate a non-linear and non-Gaussian consumption-based long-run risk asset pricing model, in which consumption volatility is modelled using an autoregressive gamma process instead of an autoregressive process that is usually adopted in the standard long-run risk model \citep*[see,][]{bansal2004risks, bansal2012empirical}. This model has also been studied by fulop2020bayesian.
The model economy consists of a representative household, a final goods producing firm, a continuum of intermediate goods producing firms, and a monetary/fiscal authority. For the purpose of self-containedness, we provide the detailed model description in what follows.
{Household.} We consider a representative household who maximizes the following expected utility
subject to the budget constraint
where $\mathbb{E}_t$ denotes conditional expectation given information up to time $t$, and $\beta \in(0,1)$ is the discount factor. The household derives utility from consumption $\mathsf{C}_t$ relative to a habit shock and real money balances $\mathsf{M}_t/\mathsf{P}_t$, with $\mathsf{P}_t$ denoting the price of the final good, and derives disutility from hours worked $\mathsf{H}_t$. $\tau$ captures the household's level of risk aversion and its inverse, $1/\tau$, is the intertemporal elasticity of substitution. $\chi_M$ is a scale factor that determines steady-state real money balances. The household receives the real wage $\mathsf{W}_t$ in exchange for labor and has access to the bond market where $\mathsf{B}_t$ nominal government bonds are traded with gross interest $\mathsf{R}_t$. Furthermore, he/she receives residual real profits $\mathsf{D}_t$ from firms and has to pay lump-sum taxes $\mathsf{T}_t$. $\mathsf{SC}_t$ is the net cash inflow from trading a full set of state-contingent securities.
{Firms.} The final goods producing firms generate aggregate output $\mathsf{Y}_t$ by combining a continuum of intermediate goods $\mathsf{Y}_t(j)$ for $j\in[0,1]$. Under the assumption of perfect competition and free entry, the demand for intermediate goods with price $\mathsf{P}_t(j)$ is given by
and the price of the final good is
where $1/\nu > 1$ represents the elasticity of demand for each intermediate good.
Intermediate good $j$ is produced by a monopolist who has the following linear production technology
where $\mathsf{N}_t(j)$ is the labor input of firm $j$ and $\mathsf{A}_t$ is an exogenous productivity process that is common to all firms and evolves according to
In the above, $\mathsf{z}_t$ captures exogenous fluctuations of the technology growth rate and $\varepsilon_{z,t}\sim\mathcal{N}(0,\sigma_z^2)$.
Firms face nominal price rigidities in terms of quadratic price adjustment costs
where $\phi$ governs the price stickiness in the economy and $\pi$ is the steady-state rate of inflation $\pi_t$, defined as $\pi_t = \mathsf{P}_t/\mathsf{P}_{t-1}$. Each firm chooses its labor input $\mathsf{N}_t(j)$ and price $\mathsf{P}_t(j)$ to maximize the present value of its future profits
where $\mathsf{Q}_{t+s|t}$ is the time $t$ value of a unit of the consumption good in period $t+s$ to the household, which is treated as exogenous by the firm.
Government policies. The government consumes a stochastic fraction of aggregate output and its spending is assumed to evolve according to
where $\mathsf{g}_t$ is an exogenous process and is assumed to follow
with $\varepsilon_{g,t}\sim\mathcal{N}(0,\sigma_g^2)$.
A central bank sets the interest rate by an interest rate feedback rule
where $\varepsilon_{R,t}\sim\mathcal{N}(0,\sigma_R^2)$ is a monetary policy shock and $\mathsf{R}_t^*$ is the nominal target rate
where $\mathsf{r}$ is the steady-state real interest rate, $\pi^*$ is the target inflation rate, which coincides in equilibrium with the steady-state inflation rate $\pi$, and $Y_t^*$ is the level of output when there are no nominal rigidities ($\phi = 0$).
The government levies a lump-sum taxes to finance any shortfalls in government revenues. Its budget constraint is given by
Summary of equilibrium conditions. We express the model in terms of detrended variables, $\mathsf{c}_t = \mathsf{C}_t/\mathsf{A}_t$ and $\mathsf{y}_t = \mathsf{Y}_t/\mathsf{A}_t$. It can be shown that if the innovations $(\varepsilon_{z,t},\varepsilon_{R,t},\varepsilon_{g,t})$ are zero at all times, the model economy has a unique steady-state, which is given as follows
Denoting $\hat{\mathsf{x}}_t = \log(\mathsf{x}_t/\mathsf{x})$ as the percentage deviation of a variable $\mathsf{x}_t$ from its steady-state $\mathsf{x}$, the model's equilibrium conditions can be summarized as follows
Model solution. The above equilibrium conditions form a non-linear rational expectations system in the state variable
in which the first four are endogenous state variables and the last three are exogenous state variables. This system has to be solved numerically before estimating the model using data; the resulting solution has the form in Equation (ref). In empirical studies, linear approximation methods sims2002solving are very popular as they result in a linear state-space representation of the model that can be easily estimated by evaluating the likelihood function using the Kalman filter an2007bayesian, herbst2016.
However, as shown in \citet*{fernandez2006ecta}, the second-order approximation errors in the solutions of dynamic economic models have first-order effects on the likelihood estimation. Furthermore, errors in the likelihood estimation are compounded with the size of the sample. Therefore, it is desirable to solve the model using more accurate methods. For the purpose of illustrating our new econometric method, after taking both accuracy and computational cost in account, we choose a second-order perturbation method with pruning \citep*{schmitt2004, schmitt2007optimal, fernandez2018res} to solve the above non-linear rational expectations system.
Measurement equations. Finally, we complete the model by relating the state variables $s_t$ to a set of observables as in Equation (ref). Following an2007bayesian and herbst2016, we assume that the following observations are available: quarter-to-quarter per capita GDP growth rates (YGR), annualized quarter-to-quarter inflation rates (INF), and annualized nominal interest rates (INT), which are in percentages and relate to the state variables as follows
where $u_t = (u_{1,t}, u_{2,t}, u_{3,t})\sim\mathcal{N}(0_3,\Sigma_u)$, which captures the measurement errors of the observables, has a zero mean vector $0_3$ and diagonal covariance matrix $\Sigma_u = \textrm{diag}(\sigma_{e,y}^2, \sigma_{e,\pi}^2, \sigma_{e,R}^2)$. The $d_{\theta}=18$ unknown parameters $\theta$ to be inferred include the structural parameters and the standard deviations of the measurement errors
where $\kappa=\tau(1-\nu)/\nu\pi^2\phi$, and $(\mathsf{r}^{(A)},\pi^{(A)},\gamma^{(Q)})$ are related to the model steady-states via
State-space model. Solving the model using the second-order perturbation method as discussed above introduces non-linearities in the resulting state-space model. The latent state variables $s_t = (x_t,z_t)\in \mathbb{S}=\mathbb{R}^{d_x}\times\mathbb{R}^{d_z}$ contain $d_x=4$ endogenous variables $x_t$ and $d_z=3$ exogenous variables $z_t$. The time evolution of the endogenous variables $(x_t)_{t=0}^T$ is given by a deterministic mapping $x_t = X_{\theta}(x_{t-1},z_t)$ for $t=1,\ldots,T$ with initialization at $x_0=X_{\theta}(0_{d_x},z_0)$, where $X_{\theta}:\mathbb{S}\rightarrow\mathbb{R}^{d_x}$ has the form
In the above, the constant is $c(\theta)\in\mathbb{R}^{d_x}$, the linear term is $L(\theta)\in\mathbb{R}^{d_x\times d}$, and the quadratic term $Q_{\theta}:\mathbb{S}\rightarrow\mathbb{R}^{d_x}$ is defined as $Q_{\theta}(s)_i=s^\top Q_i(\theta)s$ where $Q_i(\theta)\in\mathbb{R}^{d\times d}$ for $i=1,\ldots,d_x$. The exogenous variables $(z_t)_{t=0}^T$ are modelled as a vector autoregressive process
with initialization $z_0 = \Sigma(\theta)\varepsilon_0, \varepsilon_0\sim\mathcal{N}(0_{d_z},I_{d_z})$, where the matrices are $\rho(\theta),\Sigma(\theta)\in\mathbb{R}^{d_z\times d_z}$ and $I_{d_z}$ denotes the identity matrix of size $d_z$. We can then write the time evolution of the latent states $(s_t)_{t=0}^T$ as
for $t=1,\ldots,T$, with the initial state
Expressions for the model matrices $A(\theta)\in\mathbb{R}^{d\times d}$, $B(\theta)\in\mathbb{R}^{d\times d_z}$ and the function $c_{\theta}:\mathbb{S}\times\mathbb{R}^{d_z}\rightarrow\mathbb{S}$ in terms of the above quantities are given in the Appendix (ref). In the framework of Equation (ref), this corresponds to an initial distribution and a Markov transition kernel that are partially degenerate
for $t=1,\ldots,T$, where $s_0=\Phi_{\theta}^{(0)}(\varepsilon_0)$ and $s_t=\Phi_{\theta}(s_{t-1},\varepsilon_{t})$ denote the mappings in Equations (ref) and (ref). The log-linearized model can be easily obtained by simply suppressing the $c_{\theta}$ terms in Equations (ref) and (ref).
There are three observables $y_t\in\mathbb{Y}=\mathbb{R}^{d_y}$ with $d_y=3$. For each quarter $t=1,\ldots,T$, the relation between observables and latent states is described by the model
for some model vector $d(\theta)\in\mathbb{R}^{d_y}$, model matrix $E(\theta)\in\mathbb{R}^{d_y\times d}$ and covariance matrix $F(\theta)\in\mathbb{R}^{d_y\times d_y}$; see Equations (ref), (ref), and (ref).
Annealed controlled SMC. The uncontrolled SMC within AC-SMC is taken as the BPF; see Section (ref) for the corresponding choice of proposals $(q_t)_{t=0}^T$ and weight functions $(w_t)_{t=0}^T$. Owing to the partially degenerate nature of the model's initial distribution and Markov transitions in Equation (ref), here we elaborate how to specify policies. The treatment we propose is necessary for models with both endogenous and exogenous state variables and to the best of our knowledge, has not been considered in earlier works.
Instead of adopting a parameterization purely in terms of the latent variables $(s_t)_{t=0}^T$, we will also introduce dependence on the noise variables $(\varepsilon_t)_{t=0}^T$ in the function classes
for $t=1,\ldots,T$, where $Q_0(z;\beta)=z^\top A z + z^\top b + c$ is a quadratic function with coefficients $\beta=(A,b,c)\in\mathbb{R}^{d_z\times d_z}_{\mathrm{sym}}\times\mathbb{R}^{d_z}\times\mathbb{R}$ and $Q(s_{t-1},\varepsilon_t;\beta_t)$ has the form of
with coefficients $\beta=(A,b,C,D,e,f)\in\mathbb{R}^{d_z\times d_z}_{\mathrm{sym}}\times\mathbb{R}^{d_z}\times \mathbb{R}^{d_z\times d}\times \mathbb{R}_\mathrm{sym}^{d\times d}\times\mathbb{R}^d\times\mathbb{R}$. The notation $\mathbb{R}_\mathrm{sym}^{d\times d}$ refers to the space of real symmetric matrices of size $d\times d$, and $|\beta|$ denotes the Euclidean norm of $\beta$ as a vector in Euclidean space. Under this specification, we can initialize the policy in AC-SMC by having all coefficients equal to zero. Functions in Equation (ref) can be fitted using ridge regression in the logarithmic scale, and $\xi\in(0,\infty)$ controls the amount of shrinkage. Since policy approximations ultimately construct SMC proposal distributions, we take the view that careful selection of shrinkage within ADP is unnecessary. In our numerical implementation, we will fix $\xi$ as a suitably large value to safeguard against ill-conditioned Gram matrices.
The choice of function classes in Equation (ref) is well-specified when log-linearization approximations are used to solve for the equilibrium conditions; see the Appendix (ref) for the optimal policy of the linearized model. When we employ second-order approximation methods that induce non-linearities in the model, this choice can provide good policy approximations despite being misspecified. However, our numerical findings reveal that these function classes are too flexible if the coefficients are not appropriately constrained. In practice, overfitting at each step of the ADP algorithm (Algorithm (ref)) could result in numerical instabilities after several steps of iterating overfitted functions. To prevent overfitting, instead of performing regression using the variables $(s_{t-1},\varepsilon_t)$, we reduce the dimension of the feature space by working with the linearized state $\tilde{s}_t = A(\theta)s_{t-1} + B(\theta)\varepsilon_t$, and fitting a quadratic function $Q_0(\tilde{s}_t;\tilde{\beta}_t)$ with coefficients $\tilde{\beta}_t=(\tilde{A}_t,\tilde{b}_t,\tilde{c}_t)\in\mathbb{R}_\mathrm{sym}^{d\times d}\times\mathbb{R}^d\times\mathbb{R}$. The desired coefficients $\beta_t=(A_t,b_t,C_t,D_t,e_t,f_t)$ are then obtained by equating coefficients in the equality $Q_0(\tilde{s}_t;\tilde{\beta}_t)=Q(s_{t-1},\varepsilon_t;\beta_t)$, which defines the mapping $\Lambda_{\theta}$ (see, Appendix (ref) for details). We note that regularization is necessary in this setting, as the linearized states (and hence its features) are supported on a lower-dimensional subspace.
After fitting a policy $\psi=(\psi_t)_{t=0}^T$, the new proposal transitions $(q_t^{\psi})_{t=1}^T$ involve sampling the noise $\varepsilon_t$ from a distribution that depends on the previous state $s_{t-1}$, and applying the map $s_t=\Phi_{\theta}(s_{t-1},\varepsilon_{t})$ to obtain the next state (similarly for the new initial distribution $q_0^{\psi}$). As the map $\Phi_{\theta}(s_{t-1},\varepsilon_{t})$ may be the result of a complex numerical routine, the state transition density is generally intractable; hence including the noise variables in our policy parameterization is key to facilitate sampling from the new proposal transitions. We refer readers to the Appendix (ref) for the precise form of the proposal transitions and analytical expressions to evaluate the expectations of the policy appearing in the weight functions of Equation (ref) and Algorithm (ref). Lastly, as the form of the functions in Equation (ref) are closed under multiplication, policy refinements can be performed by updating the coefficients (see, Appendix (ref)).
We perform Monte Carlo simulations to examine the likelihood estimator of AC-SMC. We first consider the log-linearized model that is linear and Gaussian. Therefore, the likelihood function can be evaluated exactly using a Kalman filter. We generate a sequence of quarterly observations of YGR, INF and INT from this model for 500 time periods. The true model parameter values are the same as those used in an2007bayesian. For the purpose of comparison, we also include simulation results for the BPF with various numbers of particles. It is well-understood that the performance of the BPF can be very poor when the standard deviations of the measurement errors are small. To investigate how the efficiency of AC-SMC behaves with respect to the magnitude of the standard deviations of the measurement errors, we choose these standard deviations to be 5%, 10%, 15%, or 20% of the standard deviations of the simulated data and add measurement noise to the data accordingly. For each measurement error setting, we implement 100 independent repetitions of BPF and AC-SMC on a fixed simulated data sample and compute the sample means and variances of their log-likelihood estimates.
Table (ref) summarizes the performance of these particle filters. For AC-SMC, we set the number of particles as $N=1,024$, and for BPF, we choose the number of particles $N$ equal to 4,096, 8,192, or 16,384. It is striking that for all measurement error settings, the sample mean of the AC-SMC log-likelihood estimates matches the true log-likelihood value computed using a Kalman filter, and the sample variance is very small, particularly in the case of larger measurement errors. In contrast, these numerical results show that the BPF log-likelihood estimator is very inefficient, and especially so in the case of small measurement errors. For example, in the setting where the standard deviations of the measurement errors are 5% of their sample analogue, the true log-likelihood value is -2,091.4; the sample mean and variance of log-likelihood estimates are -2,091.4 and $2.3\times 10^{-4}$ respectively for AC-SMC, and -2,094.6 and 416.9 respectively for BPF with $N=16,384$ number of particles. In terms of computational cost, AC-SMC with $N=1,024$ number of particles is comparable to that of BPF with $N=8,192$ number of particles.
We now consider the non-linear model that is obtained by solving the DSGE model with the second-order perturbation method. Note that the likelihood function of this non-linear model is intractable. As before, we simulate data from the model and consider four measurement error settings. Table (ref) reports the sample means and variances of log-likelihood estimates from 100 independent repetitions of BPF and AC-SMC. In this simulation study, we set the number of particles in AC-SMC as $N=1,024$, and consider the number of particles $N$ equal to 16,384, 32,768, or 65,536 in BPF.
The performance of AC-SMC is also significantly better than that of BPF for this non-linear model. For example, in the challenging regime of small measurement errors, the sample variance of the log-likelihood estimates of AC-SMC is as small as 0.065, compared to 182.7 for BPF with $N=65,536$ particles. Moreover, the gains in terms of variance reduction is always of several orders of magnitude across the four measurement error settings. To take computational cost into account, we observe that the cost of AC-SMC with $N=1,024$ particles is comparable to that of BPF with a number of particles between 32,768 and 65,536. Lastly, we note that the likelihood estimate obtained by approximating the non-linear model with log-linearization and running the Kalman filter is very biased compared to the estimates from BPF and AC-SMC.
We now apply our adaptive SMC$^2$ algorithm with AC-SMC nested to estimate the above new Keynesian DSGE model using real data. While the log-linearized approximation of this model has been extensively investigated an2007bayesian, herbst2016, the corresponding non-linear model has not been fully studied yet in the literature. We estimate the resulting non-linear model from applying the second-order perturbation method with pruning based on three observables as discussed in Equations (ref), (ref), and (ref): quarterly per capita GDP growth rate (YGR), quarterly inflation (INF), and the annualized federal funds rate as a proxy of interest rate (INT). To avoid the issue of having a zero lower bound on the interest rate, we focus on the pre-crisis sample, ranging from 1983:Q1 to 2007:Q4 with a total of 100 observations. The time series of these three observables are displayed in Figure (ref).
Our adaptive SMC$^2$ algorithm requires initializing particles from the prior distribution of the model parameters. Given the complex behaviour of the likelihood function implied by the DSGE model, choosing a conjugate prior is not feasible here. We adopt a prior distribution that is component-wise independent, with marginal distributions that are normal distributions for real-valued parameters, truncated Normal distributions for positive parameters, and uniform distributions for bounded parameters. Under these choices, simulating from the prior distribution is straightforward. The hyper-parameters of the prior distribution are selected by following the existing literature in macroeconomics. The left column of Table (ref) details the exact distributional form, the support, and the hyper-parameters of the prior distribution for each model parameter.
We estimate the model in two settings: when the standard deviations of the measurement errors are fixed at 20% of the corresponding sample standard deviations, as is commonly considered in the literature, and when they are treated as free parameters to be inferred using the data. The algorithmic configuration of our adaptive SMC$^2$ involves $P=1,024$ parameter particles, $N=1,024$ state particles within AC-SMC, and an ESS threshold of $\kappa_{\mathrm{ESS}}=0.5$ to adaptively determine the inverse temperature schedule. To ensure adequate sample diversity after resampling, we apply PMMH moves until the cumulative acceptance rate reaches two. Figure (ref) displays the acceptance rate of the final PMMH move at each annealing iteration. We observe acceptance rates ranging from around 0.2 to 0.4 in both settings, revealing adequate rejuvenation of the sample diversity doucet2015efficient,sherlock2015efficiency.
The middle column of Table (ref) reports the parameter estimates of the model when standard deviations of the measurement errors are fixed. The risk-aversion coefficient has a posterior mean of 2.07 and a posterior standard deviation of 0.42. As the posterior mean of $\kappa$ is 0.76, which is relatively large, this suggests a low degree of price rigidity and a small effect of monetary policy shocks on output. The annualized steady-state growth rate of the economy is 3.2%, the annualized steady-state inflation rate is 4.8%, and the annualized steady-state nominal interest rate is 8.4%. The steady-state ratio of $\mathsf{c}/\mathsf{y}$, which is equal to $1/\mathsf{g}$, is around 0.28. The two exogenous processes, given by the government spending shock and the productivity shock, are close to unit root as the estimated persistence parameters are 0.98 and 0.99, respectively. This suggests that innovations to these processes have a long-lasting effect.
We now examine how the parameter estimates change when the standard deviations of the measurement errors are treated as free parameters. From the right column of Table (ref), the most striking change is that the posterior mean of $\kappa$ drops from 0.76 to 0.06, which implies a much larger degree of price rigidity and a stronger effect of monetary policy shocks on output. The annualized steady-state growth rate of the economy is 2.7%, the annualized steady-state inflation rate is 3.7%, the annualized steady-state nominal interest rate is 6.7%, and the steady-state ratio of $\mathsf{c}/\mathsf{y}$ is around 0.53. The government spending shock and the productivity shock are now slightly less persistent, with much smaller variations.
The estimated standard deviations of the measurement errors in the right column are quite different from the fixed values given in the middle column, and suggest that the model fits the interest rate better than the output growth rate and the inflation rate. Comparing the estimated model evidence in the two settings, reported in the last row of Table (ref), suggests that the common practice of fixing the standard deviations of the measurement errors at 20% of their sample analogue does deteriorate the model fit.
Using the output of our adaptive SMC$^2$ algorithm, we can examine the smoothing distribution of the latent states, such as the exogenous states which are usually used in policy analysis. Figure (ref) illustrates the (5, 50, 95)%-quantiles of one smoothed endogenous state, the nominal interest rate $\hat{\mathsf{R}}_{t}$, and of one smoothed exogenous state, the productivity shock $\hat{\mathsf{z}}_t$. The left and right panels correspond to fixing the standard deviations of the measurement errors and treating them as free parameters, respectively. The differences between the smoothing distributions of these state variables are apparent.
The literature on estimation and empirical applications of dynamic structural macrofinance models has mostly been relying on approximations using log-linearization to cast these models into linear Gaussian state-space models. It is now known that second-order approximation errors in the linear solutions can have first-order effects on likelihood estimation, and that errors in likelihood estimation can be compounded with the sample size. However, when models are solved using more accurate numerical methods such as high-order perturbation methods or projection methods, the solved models become non-linear and/or non-Gaussian state-space models with high-dimensional and complex structures.
In this paper, we propose a novel methodology to efficiently estimate the likelihood function of such state-space models by building on state-of-the-art SMC methods. While particle filters have been applied to estimate the likelihood of dynamic economic models, the successful application of this approach remains challenging due to the large variance of the resulting likelihood estimator. We develop an annealed controlled SMC method that delivers numerically stable and low variance estimators of the likelihood function, by adopting an annealing procedure to gradually introduce information from observations and construct globally optimal proposal distributions using approximate dynamic programming schemes. We provide a theoretical analysis to characterize the quality of our proposal distributions, which elucidates various properties of annealed controlled SMC.
To perform parameter inference, we develop a new adaptive SMC$^2$ algorithm that employs likelihood estimators from annealed controlled SMC. Under suitable assumptions, we establish law of large numbers and central limit theorems for our estimators of posterior expectations and the model evidence. We illustrate the strengths of our adaptive SMC$^2$ algorithm by estimating two popular macrofinance models: a prototypical new Keynesian dynamic stochastic general equilibrium model and a non-linear non-Gaussian consumption-based long-run risk asset pricing model.
We believe that the methods developed in this paper are readily applicable to DSGE models exhibiting strong non-linearities for instance due to the presence of effective lower bounds on interest rates targeted by monetary policy. Further as researchers move to the estimation of larger scale non-linear DSGE models with richer state spaces, the curse of dimensionality of particle filters (see e.g. bengtsson2008curse) will call for the use of more efficient approaches, such as the one developed here.
This work was funded by CY Initiative of Excellence (grant "Investissements d'Avenir" ANR-16-IDEX-0008).