# Efficient Likelihood-based Estimation via Annealing for Dynamic Structural
Macrofinance Models

Andras Fulop ESSEC Business School, Paris-Singapore. Email<EMAIL_ADDRESS>Jeremy Heng ESSEC Business School, Paris-Singapore. Email<EMAIL_ADDRESS>Junye Li Fudan University. Email<EMAIL_ADDRESS>
(First version: January 2021; This version: December 2021)

###### Abstract

Most solved dynamic structural macrofinance models are non-linear and/or non-
Gaussian state-space models with high-dimensional and complex structures. We
propose an annealed controlled sequential Monte Carlo method that delivers
numerically stable and low variance estimators of the likelihood function. The
method relies on an annealing procedure to gradually introduce information
from observations and constructs globally optimal proposal distributions by
solving associated optimal control problems that yield zero variance
likelihood estimators. To perform parameter inference, we develop a new
adaptive SMC2 algorithm that employs likelihood estimators from annealed
controlled sequential Monte Carlo. We provide a theoretical stability analysis
that elucidates the advantages of our methodology and asymptotic results
concerning the consistency and convergence rates of our SMC2 estimators. We
illustrate the strengths of our proposed methodology by estimating two popular
macrofinance models: a non-linear new Keynesian dynamic stochastic general
equilibrium model and a non-linear non-Gaussian consumption-based long-run
risk model.

Keywords: Sequential Monte Carlo, Particle Filters, Approximate Dynamic
Programming, Annealing, SMC2, DSGE, Long-Run Risk.

##  Introduction

In macroeconomics and finance, dynamic structural models provide a convenient
theoretical framework to explain economic fluctuations, government policies,
and asset prices movements.111For example, dynamic stochastic general
equilibrium (DSGE) models, spanning from real business cycle models (Kydland
and Prescott, 1982) to new Keynesian models (see, e.g., Christiano,
Eichenbaum, and Evans, 2005; Smets and Wouters, 2003), 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 (Lucas, 1978). 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.222The 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, Herbst and Schorfheide (2016) and
Fernandez-Villaverde, Rubio-Ramirez, and Schorfheide (2016). 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, Bansal, Gallant, and Tauchen (2007), Bansal, Kiku,
and Yaron (2016a), and Gallant, Jahan-Parvar, and Liu (2019).

It is now known that second-order approximation errors in linear solutions of
dynamic economic models have first-order effects on likelihood estimation
(Fernandez-Villaverde, Rubio-Ramirez, and Santos, 2006; Fernandez-Villaverde
and Rubio-Ramirez, 2007), and that errors in the likelihood estimation are
compounded with the sample size. In recent work, Pohl, Schmedders, and Wilms
(2018) 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 (Schmitt-Grohe and
Uribe, 2004) or projection methods (Judd, 1992), 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 (see, e.g.,
Fernandez-Villaverde and Rubio-Ramirez, 2007; DeJong et al., 2013; Herbst and
Schorfheide, 2019; Fulop et al., 2021), 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 Heng et al. (2020) 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 Heng et al. (2020) 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 (Herbst and Schorfheide, 2019; Godsill and Clapp,
2001). In an offline context where globally optimal proposals are desired, the
use of annealing to alleviate numerical instabilities was mentioned in earlier
work by Scharth and Kohn (2016), 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 SMC2 algorithm that recursively approximates a
sequence of annealed posterior distributions of the latent states and model
parameters. Compared to existing annealed SMC2 routines such as Duan and Fulop
(2015) and Svensson, Schön, and Lindsten (2018), 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 (Andrieu, Doucet, and
Holenstein, 2010) that sample from the posterior distribution using a single
Markov chain, our algorithm shares the many advantages of SMC samplers (Del
Moral et al., 2006; Dai et al., 2020). 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 SMC2 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 Woodford (2003), An and
Schorfheide (2007) and Herbst and Schorfheide (2016); the second is a
consumption-based long-run risk asset pricing model, in which consumption
volatility is modelled as an autoregressive gamma process (Gourieroux and
Jasiak, 2006) instead of an autoregressive process that is typically adopted
in the standard long-run risk model (see, e.g., Bansal and Yaron, 2004;
Bansal, Kiku, and Yaron, 2012a). 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 (Sims, 2002) or the
second-order perturbation method with pruning (Schmitt-Grohe and Uribe, 2004,
2007; Andreasen, Fernandez-Villaverde, and Rubio-Ramirez, 2018). 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 (Gordon et al., 1993), 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 An and Schorfheide (2007)
and empirically estimated in Herbst and Schorfheide (2016), the corresponding
non-linear model has not been fully studied yet in the literature.333In a
recent paper by Aruoba et al. (2021), 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 2 describes the
models of interest and its state-space representation. Section 3 introduces
our annealed controlled SMC method to efficiently estimate the likelihood of
these models. Section 4 details our adaptive SMC2 algorithm for parameter
inference. Section 5 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 6 concludes the paper. Supplementary details and proofs of all
theoretical results are provided in the Appendix.

##  State-space representation

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

$\displaystyle s_{0}=\Phi_{\theta}^{(0)}(\varepsilon_{0}),\quad
s_{t}=\Phi_{\theta}(s_{t-1},\varepsilon_{t}),\quad
y_{t}=\Psi_{\theta}(s_{t-1},s_{t},u_{t}),$ (1)

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

$s_{0}\sim\mu_{\theta},\quad s_{t}\sim f_{\theta}(\cdot|s_{t-1}),\quad
t=1,\ldots,T,$ (2)

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

$y_{t}\sim g_{\theta}(\cdot|s_{t-1},s_{t}),\quad t=1,\ldots,T,$ (3)

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 (3) is convenient when modeling asset
returns data; longer time dependencies can also be accommodated in our
proposed methodology.

##  Likelihood estimation

Given an observation sequence $y_{1:T}=(y_{t})_{t=1}^{T}\in\mathbb{Y}^{T}$,
the complete likelihood is

$p(ds_{0:T},y_{1:T}|\theta)=p(ds_{0:T}|\theta)p(y_{1:T}|s_{0:T},\theta)=\mu_{\theta}(ds_{0})\prod_{t=1}^{T}f_{\theta}(ds_{t}|s_{t-1})g_{\theta}(y_{t}|s_{t-1},s_{t}).$
(4)

To perform parameter inference, we have to integrate out the latent process to
compute the likelihood function

$p(y_{1:T}|\theta)=\int_{\mathbb{S}^{T+1}}p(ds_{0:T},y_{1:T}|\theta),$ (5)

and to estimate the latent states, we need to characterize the smoothing
distribution

$p(ds_{0:T}|y_{1:T},\theta)=\frac{p(ds_{0:T},y_{1:T}|\theta)}{p(y_{1:T}|\theta)}.$
(6)

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

$\displaystyle p(ds_{0:T},y_{1:T}|\theta,\lambda)$ $\displaystyle=$
$\displaystyle p(ds_{0:T}|\theta)p(y_{1:T}|s_{0:T},\theta,\lambda)$ (7)
$\displaystyle=$
$\displaystyle\mu_{\theta}(ds_{0})\prod_{t=1}^{T}f_{\theta}(ds_{t}|s_{t-1})g_{\theta}(y_{t}|s_{t-1},s_{t})^{\lambda},$

and the corresponding likelihood and smoothing distribution as

$p(y_{1:T}|\theta,\lambda)=\int_{\mathbb{S}^{T+1}}p(ds_{0:T},y_{1:T}|\theta,\lambda),\quad
p(ds_{0:T}|y_{1:T},\theta,\lambda)=\frac{p(ds_{0:T},y_{1:T}|\theta,\lambda)}{p(y_{1:T}|\theta,\lambda)}.$
(8)

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 (6). In the
context of state-space models, similar motivations can be found in Godsill and
Clapp (2001), Svensson, Schön, and Lindsten (2018), and Herbst and Schorfheide
(2019). We note that the likelihood $p(y_{1:T}|\theta,\lambda)$ in Equation
(8) is different but related to the tempered likelihood
$p(y_{1:T}|\theta)^{\lambda}$, considered in Duan and Fulop (2015), via
Jensen’s inequality.

###  Sequential Monte Carlo

For most models of practical interest, the quantities in Equation (8) 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}$ (Doucet, De Freitas, and Gordon, 2001; Chopin and
Papaspiliopoulos, 2020). 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}^{N}W_{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

$q(ds_{0:T}|\theta,\lambda)=q_{0}(ds_{0}|\theta,\lambda)\prod_{t=1}^{T}q_{t}(ds_{t}|s_{t-1},\theta,\lambda)$
(9)

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

$w_{0}(s_{0};\theta,\lambda)\prod_{t=1}^{T}w_{t}(s_{t-1},s_{t};\theta,\lambda)=\frac{p(s_{0:T},y_{1:T}|\theta,\lambda)}{q(s_{0:T}|\theta,\lambda)}.$
(10)

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 (Jacob,
Murray, and Rubenthaler, 2015). 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

$\hat{p}(y_{1:T}|\theta,\lambda)=\left\\{\frac{1}{N}\sum_{n=1}^{N}w_{0}(s_{0}^{(n)};\theta,\lambda)\right\\}\left\\{\prod_{t=1}^{T}\frac{1}{N}\sum_{n=1}^{N}w_{t}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)};\theta,\lambda)\right\\},$
(11)

and a weighted particle approximation of the corresponding smoothing
distribution is given by

$\hat{p}(ds_{0:T}|y_{1:T},\theta,\lambda)=\sum_{n=1}^{N}W_{T}^{(n)}\delta_{s_{0:T}^{(n)}}(ds_{0:T}).$
(12)

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 (12), 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}^{N}W_{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 (Del Moral, 2004; Chopin, 2004). However, a successful
implementation of SMC in practice crucially relies on the choice of proposal
distributions in Equation (9). 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 Gordon, Salmond, and
Smith (1993) corresponds to having the law of the latent process in Equation
(2) 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 (10). 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 Pitt and Shephard (1999) and Carpenter, Clifford, and
Fearnhead (1999) 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
(see, e.g., Doucet et al., 2000; Finke et al., 2020).

###  Controlled sequential Monte Carlo

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 Heng et al.
(2020) that in turn builds upon the works by Richard and Zhang (2007); Scharth
and Kohn (2016) and Guarniero, Johansen, and Lee (2017).

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 (10). 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}^{T}p(ds_{t}|s_{t-1},y_{t:T},\theta,\lambda)$
as

$\displaystyle
p(ds_{0}|y_{1:T},\theta,\lambda)=\frac{q_{0}(ds_{0})\psi_{0}^{*}(s_{0})}{q_{0}(\psi_{0}^{*})},\quad
p(ds_{t}|s_{t-1},y_{t:T},\theta,\lambda)=\frac{q_{t}(ds_{t}|s_{t-1})\psi_{t}^{*}(s_{t-1},s_{t})}{q_{t}(\psi_{t}^{*}|s_{t-1})}$
(13)

for $t=1,\ldots,T$, where the sequence of functions
$\psi^{*}=(\psi_{t}^{*})_{t=0}^{T}$ are defined by the backward recursion

$\displaystyle\psi_{T}^{*}(s_{T-1},s_{T})~{}$
$\displaystyle=~{}w_{T}(s_{T-1},s_{T}),$
$\displaystyle\psi_{t}^{*}(s_{t-1},s_{t})~{}$
$\displaystyle=~{}w_{t}(s_{t-1},s_{t})q_{t+1}(\psi_{t+1}^{*}|s_{t}),\quad
t=T-1,\ldots,1,$ (14) $\displaystyle\psi_{0}^{*}(s_{0})~{}$
$\displaystyle=~{}w_{0}(s_{0})q_{1}(\psi_{1}^{*}|s_{0}).$

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

$\psi_{0}^{*}(s_{0})=p(y_{1:T}|s_{0},\theta,\lambda),\quad\psi_{t}^{*}(s_{t-1},s_{t})=p(y_{t:T}|s_{t-1},s_{t},\theta,\lambda),\quad
t=1,\ldots,T,$ (15)

it is sometimes referred to as the backward information filter (Briers,
Doucet, and Maskell, 2010).

By obtaining an approximation $\psi=(\psi_{t})_{t=0}^{T}$ of the backward
recursion in Equation (3.2) using an approximate dynamic programming method
that will be discussed in Section 3.3, we can construct a new proposal
distribution

$\displaystyle
q^{\psi}(ds_{0:T})=q_{0}^{\psi}(ds_{0})\prod_{t=1}^{T}q_{t}^{\psi}(ds_{t}|s_{t-1})$
(16)

by mimicking Equation (13), i.e. define

$\displaystyle
q_{0}^{\psi}(ds_{0})=\frac{q_{0}(ds_{0})\psi_{0}(s_{0})}{q_{0}(\psi_{0})},\quad
q_{t}^{\psi}(ds_{t}|s_{t-1})=\frac{q_{t}(ds_{t}|s_{t-1})\psi_{t}(s_{t-1},s_{t})}{q_{t}(\psi_{t}|s_{t-1})},\quad
t=1,\ldots,T.$ (17)

Following the terminology in Heng et al. (2020), 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 5 and assume $\psi$ is such
that the proposals $(q_{t}^{\psi})_{t=0}^{T}$ in Equation (17) 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

$\displaystyle w_{0}^{\psi}(s_{0})~{}$
$\displaystyle=~{}\frac{q_{0}(\psi_{0})q_{1}(\psi_{1}|s_{0})}{\psi_{0}(s_{0})},$
$\displaystyle w_{t}^{\psi}(s_{t-1},s_{t})~{}$
$\displaystyle=~{}\frac{w_{t}(s_{t-1},s_{t})q_{t+1}(\psi_{t+1}|s_{t})}{\psi_{t}(s_{t-1},s_{t})},\quad
t=1,\ldots,T-1,$ (18) $\displaystyle w_{T}^{\psi}(s_{T-1},s_{T})~{}$
$\displaystyle=~{}\frac{w_{T}(s_{T-1},s_{T})}{\psi_{T}(s_{T-1},s_{T})},$

which satisfy
$w_{0}^{\psi}(s_{0})\prod_{t=1}^{T}w_{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 1. 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 1, 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 (12) (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
3.5. We note that the choice $\psi=\psi^{*}$ is optimal as Algorithm 1 would
yield a zero variance likelihood estimator, for any number of particles $N$,
and an exact trajectory from the smoothing distribution.

Algorithm 1 Controlled sequential Monte Carlo at inverse temperature
$\lambda\in[0,1]$

Input: number of particles $N$ and policy $\psi=(\psi_{t})_{t=0}^{T}$.

(1) For time $t=0$ and particle $n=1,\ldots,N$.

(1a) Sample state $s_{0}^{(n)}\sim q_{0}^{\psi}$.

(1b) Compute normalized weights
$W_{0}^{(n)}=w_{0}^{\psi}(s_{0}^{(n)})/\sum_{m=1}^{N}w_{0}^{\psi}(s_{0}^{(m)})$.

(2) For time $t=1,\ldots,T$ and particle $n=1,\ldots,N$.

(2a) Sample ancestor $a_{t-1}^{(n)}\sim
r(\cdot|W_{t-1}^{(1)},\ldots,W_{t-1}^{(N)})$.

(2b) Sample state $s_{t}^{(n)}\sim
q_{t}^{\psi}(\cdot|s_{t-1}^{(a_{t-1}^{(n)})})$.

(2c) Compute normalized weights
$W_{t}^{(n)}=w_{t}^{\psi}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})/\sum_{m=1}^{N}w_{t}^{\psi}(s_{t-1}^{(a_{t-1}^{(m)})},s_{t}^{(m)})$.

(3) Compute likelihood estimator
$\hat{p}(y_{1:T}|\theta,\lambda)=\\{\frac{1}{N}\sum_{n=1}^{N}w_{0}^{\psi}(s_{0}^{(n)})\\}\\{\prod_{t=1}^{T}\frac{1}{N}\sum_{n=1}^{N}w_{t}^{\psi}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})\\}$.

(4) Sample an ancestor $l_{T}\sim r(\cdot|W_{T}^{(1)},\ldots,W_{T}^{(N)})$ and
set ancestral lineage as $l_{t}=a_{t}^{(l_{t+1})}$ for $t=T-1,\ldots,0$.

Output: states $(s_{t}^{(n)})_{t=0,n=1}^{T,N}$, ancestors
$(a_{t}^{(n)})_{t=0,n=1}^{T-1,N}$, likelihood estimator
$\hat{p}(y_{1:T}|\theta,\lambda)$ and trajectory
$(s_{t}^{(l_{t})})_{t=0}^{T}$.

###  Policy refinement with approximate dynamic programming

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 (3.2). If we define $\phi^{*}=(\phi_{t}^{*})_{t=0}^{T}$
using the backward recursion

$\displaystyle\phi_{T}^{*}(s_{T-1},s_{T})$ $\displaystyle=$ $\displaystyle
w_{T}^{\psi}(s_{T-1},s_{T}),$ $\displaystyle\phi_{t}^{*}(s_{t-1},s_{t})$
$\displaystyle=$ $\displaystyle
w_{t}^{\psi}(s_{t-1},s_{t})q_{t+1}^{\psi}(\phi_{t+1}^{*}|s_{t}),\quad
t=T-1,\ldots,1,$ (19) $\displaystyle\phi_{0}^{*}(s_{0})$ $\displaystyle=$
$\displaystyle w_{0}^{\psi}(s_{0})q_{1}^{\psi}(\phi_{1}^{*}|s_{0}),$

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
(see, Proposition 1, Heng et al., 2020). 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 (3.3) is the corresponding dynamic programming recursion. As noted by
Heng et al. (2020), 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 (3.3).

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

$\displaystyle\phi_{T}=\arg\min_{f\in\mathbb{F}_{T}}\sum_{n=1}^{N}\left(\log
f(s_{T-1}^{(a_{T-1}^{(n)})},s_{T}^{(n)})-\log
w_{T}^{\psi}(s_{T-1}^{(a_{T-1}^{(n)})},s_{T}^{(n)})\right)^{2},$ (20)

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

$\displaystyle\phi_{T-1}=\arg\min_{f\in\mathbb{F}_{T-1}}\sum_{n=1}^{N}\left(\log
f(s_{T-2}^{(a_{T-2}^{(n)})},s_{T-1}^{(n)})-\log\varphi_{T-1}(s_{T-2}^{(a_{T-2}^{(n)})},s_{T-1}^{(n)})\right)^{2},$
(21)

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 2, 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 2 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 (3.2), 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^{*}$.

Algorithm 2 Approximate dynamic programming at inverse temperature
$\lambda\in[0,1]$

Input: current policy $\psi=(\psi_{t})_{t=0}^{T}$ and output of controlled SMC
(Algorithm 1).

(1) For time $t=T,\ldots,1$,

(1a) If $t=T$, for particle $n=1,\ldots,N$, set
$\varphi_{T}(s_{T-1}^{(a_{T-1}^{(n)})},s_{T}^{(n)})=w_{T}^{\psi}(s_{T-1}^{(a_{T-1}^{(n)})},s_{T}^{(n)})=\frac{w_{T}(s_{T-1}^{(a_{T-1}^{(n)})},s_{T}^{(n)})}{\psi_{T}(s_{T-1}^{(a_{T-1}^{(n)})},s_{T}^{(n)})}$.

(1b) If $t<T$, for particle $n=1,\ldots,N$, set

$\varphi_{t}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})=w_{t}^{\psi}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})q_{t+1}^{\psi}(\phi_{t+1}|s_{t}^{(n)})=\frac{w_{t}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})q_{t+1}(\psi_{t+1}\cdot\phi_{t+1}|s_{t}^{(n)})}{\psi_{t}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})}.$

(1c) Fit the function
$\phi_{t}=\arg\min_{f\in\mathbb{F}_{t}}\sum_{n=1}^{N}\left(\log
f(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})-\log\varphi_{t}(s_{t-1}^{(a_{t-1}^{(n)})},s_{t}^{(n)})\right)^{2}.$

(2) For time $t=0$.

(2a) For particle $n=1,\ldots,N$, set
$\varphi_{0}(s_{0}^{(n)})=w_{0}^{\psi}(s_{0}^{(n)})q_{1}^{\psi}(\phi_{1}|s_{0}^{(n)})=q_{0}(\psi_{0})q_{1}(\psi_{1}\cdot\phi_{1}|s_{0}^{(n)})/\psi_{0}(s_{0}^{(n)}).$

(2b) Fit the function
$\phi_{0}=\arg\min_{f\in\mathbb{F}_{0}}\sum_{n=1}^{N}\left(\log
f(s_{0}^{(n)})-\log\varphi_{0}(s_{0}^{(n)})\right)^{2}.$

Output: refined policy $\psi\cdot\phi=(\psi_{t}\cdot\phi_{t})_{t=0}^{T}$.

We can then run a controlled SMC (Algorithm 1) 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 2) 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 Heng et
al. (2020). 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 (6). 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.

###  Annealed controlled sequential Monte Carlo

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 (5) and (6), approximating the
intermediate quantities in Equation (8) for $0\leq\lambda<1$ is also useful
when we consider parameter inference in Section 4.

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 (2) 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 2 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 1) at inverse temperature $\lambda_{i}$
with new proposals defined by $\psi^{(i)}$. We provide an algorithmic summary
of the above methodology in Algorithm 3, which we refer to as annealed
controlled SMC (AC-SMC).

Algorithm 3 Annealed controlled sequential Monte Carlo

Input: number of particles $N$ and inverse temperature schedule
$(\lambda_{i})_{i=0}^{I}$.

(1) For iteration $i=0$.

(1a) Set policy $\psi^{(0)}$ as constant one functions.

(1b) Run uncontrolled SMC at inverse temperature $\lambda_{0}$ with proposals
$(q_{t})_{t=0}^{T}$ and weights $(w_{t})_{t=0}^{T}$.

(2) For iteration $i=1,\ldots,I$.

(2a) Run ADP (Algorithm 2) at inverse temperature $\lambda_{i}$ with
$\psi^{(i-1)}$ as current policy and previous SMC output to obtain refined
policy $\psi^{(i)}$.

(2b) Run controlled SMC (Algorithm 1) at inverse temperature $\lambda_{i}$
with policy $\psi^{(i)}$.

Output: policy $\psi^{(I)}$ and controlled SMC output at inverse temperature
$\lambda_{I}$.

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 (3.3) and can be rewritten as

$\displaystyle\phi_{T}^{*}(s_{T-1},s_{T};\theta,\lambda_{i})~{}$
$\displaystyle=w_{T}^{\psi^{(i-1)}}(s_{T-1},s_{T};\theta,\lambda_{i-1})\frac{w_{T}(s_{T-1},s_{T};\theta,\lambda_{i})}{w_{T}(s_{T-1},s_{T};\theta,\lambda_{i-1})},$
(22) $\displaystyle\phi_{t}^{*}(s_{t-1},s_{t};\theta,\lambda_{i})~{}$
$\displaystyle=~{}w_{t}^{\psi^{(i-1)}}(s_{t-1},s_{t};\theta,\lambda_{i-1})\frac{w_{t}(s_{t-1},s_{t};\theta,\lambda_{i})}{w_{t}(s_{t-1},s_{t};\theta,\lambda_{i-1})}q_{t+1}^{\psi^{(i-1)}}(\phi_{t+1}^{*}|s_{t},\theta),$
$\displaystyle\phi_{0}^{*}(s_{0};\theta,\lambda_{i})~{}$
$\displaystyle=~{}w_{0}^{\psi^{(i-1)}}(s_{0};\theta)q_{1}^{\psi^{(i-1)}}(\phi_{1}^{*}|s_{0},\theta),$

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 Heng et al. (2020), 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 Scharth and
Kohn (2016), 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 3), 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 (Kong, Liu, and Wong, 1994), defined at each time $t=0,1,\ldots,T$
as $1/\sum_{n=1}^{N}(W_{t}^{(n)})^{2}$. Using the output of Algorithm 3, 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})$.

###  Analysis of annealed controlled sequential Monte Carlo

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 2) at inverse temperature
$\lambda\in[0,1]$. To measure the effect policy refinement has on controlled
SMC (Algorithm 1), we consider the Kullback–Leibler (KL) divergence from the
proposal distribution $q^{\psi\cdot\phi}(ds_{0:T})$, defined in Equations (16)
and (17), to the smoothing distribution $p(ds_{0:T}|y_{1:T},\theta,\lambda)$
in Equation (8), defined as

$\displaystyle\mathrm{KL}\left(p(ds_{0:T}|y_{1:T},\theta,\lambda)~{}|~{}q^{\psi\cdot\phi}(ds_{0:T})\right)=\int_{\mathbb{S}^{T+1}}\log\left(\frac{p(s_{0:T}|y_{1:T},\theta,\lambda)}{q^{\psi\cdot\phi}(s_{0:T})}\right){p(s_{0:T}|y_{1:T},\theta,\lambda)}ds_{0:T}.$
(23)

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
(Chatterjee and Diaconis, 2018), 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 A.

######  Proposition 1.

For any current policy $\psi=(\psi_{t})_{t=0}^{T}$ and an approximation
$\phi=(\phi_{t})_{t=0}$ of $\phi^{*}=(\phi_{t}^{*})_{t=0}$, the optimal
refinement of $\psi$, the KL divergence from $q^{\psi\cdot\phi}(ds_{0:T})$ to
$p(ds_{0:T}|y_{1:T},\theta,\lambda)$ satisfies

$\displaystyle\mathrm{KL}\left(p(ds_{0:T}|y_{1:T},\theta,\lambda)~{}|~{}q^{\psi\cdot\phi}(ds_{0:T})\right)\leq\sum_{t=0}^{T}\xi_{t}^{*}(\log(\phi_{t}^{*}/\phi_{t}))+\xi_{t}^{\psi\cdot\phi}(\log(\phi_{t}/\phi_{t}^{*})).$
(24)

Next, we will introduce some assumptions which are needed to derive upper
bounds of the terms in Equation (24). 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

$\displaystyle\nu_{0}^{\psi,N}(ds_{0})=\frac{1}{N}\sum_{n=1}^{N}\delta_{s_{0}^{(n)}}(ds_{0}),\quad\nu_{t}^{\psi,N}(ds_{t-1},ds_{t})=\frac{1}{N}\sum_{n=1}^{N}\delta_{(s_{t-1}^{a_{t-1}^{(n)}},s_{t}^{(n)})}(ds_{t-1},ds_{t}),$
(25)

for $t=1,\ldots,T$, that are used in the least squares approximations in
Equations (20) and (21). 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

$\displaystyle\nu_{0}^{\psi}(ds_{0})~{}$
$\displaystyle=~{}q_{0}^{\psi}(ds_{0}),\quad\nu_{1}^{\psi}(ds_{0},ds_{1})~{}=~{}\frac{\nu_{0}^{\psi}(ds_{0})w_{0}^{\psi}(s_{0})}{\nu_{0}^{\psi}(w_{0}^{\psi})}q_{1}^{\psi}(ds_{1}|s_{0}),$
(26) $\displaystyle\nu_{t}^{\psi}(ds_{t-1},ds_{t})~{}$
$\displaystyle=~{}\frac{\int_{\mathbb{S}}\nu_{t-1}^{\psi}(ds_{t-2},ds_{t-1})w_{t-1}^{\psi}(s_{t-2},s_{t-1})}{\nu_{t-1}^{\psi}(w_{t-1}^{\psi})}q_{t}^{\psi}(ds_{t}|s_{t-1}),\quad
t=2,\ldots,T.$

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})$ space444With 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 2) in
the infinite particle regime, we define $L^{2}$-projections under the function
class $\mathbb{F}_{t}$ and the distribution $\nu_{t}^{\psi}$

$\displaystyle P_{t}^{\psi}\varphi=\arg\min_{f\in\mathbb{F}_{t}}\|\log
f-\log\varphi\|_{L^{2}(\nu_{t}^{\psi})}^{2},$ (27)

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

The function classes $(\mathbb{F}_{t})_{t=0}^{T}$ satisfy:

  1. (i)

$\log\mathbb{F}_{t}$ is a closed linear subspace of $L^{2}(\nu_{t}^{\psi})$
for $t=0,\ldots,T$;

  2. (ii)

$\sup_{f\in\mathbb{G}_{t}^{\psi}}\|\log P_{t}^{\psi}f-\log
f\|_{L^{2}(\nu_{t}^{\psi})}\leq e_{t}^{\psi}<\infty$ for $t=0,\ldots,T$, where
$\mathbb{G}_{t}^{\psi}=\\{w_{t}^{\psi}q_{t+1}^{\psi}(\varphi):\varphi\in\mathbb{F}_{t+1}\\}$
for $t=0,\ldots,T-1$ and $\mathbb{G}_{T}^{\psi}=\\{w_{T}^{\psi}\\}$.

Assumption 1$(i)$ ensures the existence of a unique projection in Equation
(27), 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
1$(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 2.

There exist positive constants $(C_{t})_{t=0}^{T}$ and $(M_{t})_{t=0}^{T}$
such that:

  1. (i)

$\xi_{t}^{*}\ll\nu_{t}^{\psi}$ with density satisfying
$\xi_{t}^{*}/\nu_{t}^{\psi}\leq C_{t}$ for $t=0,\ldots,T$;

  2. (ii)

$\xi_{t}^{\psi\cdot\phi}\ll\xi_{t}^{*}$ with density satisfying
$\xi_{t}^{\psi\cdot\phi}/\xi_{t}^{*}\leq M_{t}$ for $t=0,\ldots,T$.

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

######  Theorem 1.

Under Assumptions 1 and 2,
$\varepsilon_{t}^{*}=\xi_{t}^{*}(\log(\phi_{t}^{*}/\phi_{t}))$ and
$\varepsilon_{t}^{\psi\cdot\phi}=\xi_{t}^{\psi\cdot\phi}(\log(\phi_{t}/\phi_{t}^{*}))$
for $t=0,\ldots,T$ satisfy the backward recursions

$\displaystyle\varepsilon_{t}^{*}\leq\varepsilon_{t+1}^{*}+C_{t}e_{t}^{\psi},\quad\varepsilon_{t}^{\psi\cdot\phi}\leq
M_{t}\varepsilon_{t+1}^{\psi\cdot\phi}+C_{t}M_{t}e_{t}^{\psi},\quad
t=0,\ldots,T-1,$ (28)

with $\varepsilon_{T}^{*}\leq C_{T}e_{T}^{\psi}$ and
$\varepsilon_{t}^{\psi\cdot\phi}\leq C_{T}M_{T}e_{T}^{\psi}$ at the terminal
time $T$.

Equation (28) 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 1 and
Proposition 1 would imply good performance of AC-SMC.

##  Parameter inference

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

$p(d\theta,ds_{0:T}|y_{1:T})=p(d\theta|y_{1:T})p(ds_{0:T}|y_{1:T},\theta)=\frac{p(d\theta)p(y_{1:T}|\theta)}{p(y_{1:T})}p(ds_{0:T}|y_{1:T},\theta),$
(29)

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 (Del Moral et al., 2006) to
sequentially approximate distributions and their normalizing constants along
the path

$\displaystyle p(d\theta,ds_{0:T}|y_{1:T},\lambda)$ $\displaystyle=$
$\displaystyle p(d\theta|y_{1:T},\lambda)p(ds_{0:T}|y_{1:T},\theta,\lambda)$
(30) $\displaystyle=$
$\displaystyle\frac{p(d\theta)p(y_{1:T}|\theta,\lambda)}{p(y_{1:T}|\lambda)}p(ds_{0:T}|y_{1:T},\theta,\lambda)$

for $\lambda\in[0,1]$, in which a nested AC-SMC method (Algorithm 3) 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 (30) 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 (29), 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 SMC2 methods developed by Fulop and Li
(2013), Chopin, Jacob, and Papaspiliopoulos (2013), and Duan and Fulop (2015),
but differs in important aspects. In Section 4.1, we provide a description of
our adaptive SMC2 algorithm and discuss how it differs from and relates to the
existing literature. Theoretical justifications of the algorithm are presented
in Section 4.2.

###  Adaptive SMC2

Our proposed SMC2 methodology, detailed in Algorithm 4, builds on AC-SMC in
Algorithm 3 and a conditional implementation of controlled SMC (Algorithm 5 in
Appendix B).

Algorithm 4 Adaptive sequential Monte Carlo2

Input: number of parameter particles $P$, state particles $N$ and ESS
threshold $\kappa_{\mathrm{ESS}}$.

(1) For iteration $i=0$.

(1a) Set inverse temperature as $\lambda_{0}=0$.

(1b) For particle $p=1,\ldots,P$, sample a parameter $\theta^{(p)}\sim
p(d\theta)$ from the prior and a trajectory
$s_{0:T}^{(p)}=(s_{t}^{(p)})_{t=0}^{T}\sim p(ds_{0:T}|\theta^{(p)})$ from the
latent process.

(2) For iteration $i=1,2,\ldots$, if
$\mathrm{ESS}\left(\Omega_{i}^{(1)}(1),\ldots,\Omega_{i}^{(P)}(1)\right)\geq\kappa_{\mathrm{ESS}}P$,
set next inverse temperature as $\lambda_{i}=1$, else determine $\lambda_{i}$
as the value of $\lambda\in(\lambda_{i-1},1)$ that solves
$\mathrm{ESS}\left(\Omega_{i}^{(1)}(\lambda),\ldots,\Omega_{i}^{(P)}(\lambda)\right)=\kappa_{\mathrm{ESS}}P$
using bisection method.

(3) For particle $p=1,\ldots,P$.

(3a) Compute unnormalized weights
$\omega_{i}^{(p)}=\omega_{i}^{(p)}(\lambda_{i})$ and normalized weights
$\Omega_{i}^{(p)}=\Omega_{i}^{(p)}(\lambda_{i})$.

(3b) Sample ancestor $\alpha_{i}^{(p)}\sim
r(\cdot|\Omega_{i}^{(1)},\ldots,\Omega_{i}^{(P)})$, and update parameter as
$\theta^{(p)}=\theta^{(\alpha_{i}^{(p)})}$ and trajectory as
$s_{0:T}^{(p)}=s_{0:T}^{(\alpha_{i}^{(p)})}$.

(3c) Learn policy $\psi^{(p)}$ for parameter $\theta^{(p)}$ at inverse
temperature $\lambda_{i}$.

(3d) Run conditional SMC (Algorithm 5 in the Appendix) with policy
$\psi^{(p)}$ and reference trajectory $s_{0:T}^{(p)}$ at parameter
$\theta^{(p)}$ and inverse temperature $\lambda_{i}$ to update likelihood
estimator $\hat{p}(y_{1:T}|\theta^{(p)},\lambda_{i})$ and trajectory
$s_{0:T}^{(p)}$.

(4) Construct proposal transition kernel $h_{i}$ based on current set of
parameters $(\theta^{(p)})_{p=1}^{P}$. For $k=1,\ldots,K$ PMMH moves.

(4a) Sample a parameter proposal $\theta^{(*)}\sim h_{i}(\cdot|\theta^{(p)})$.

(4b) Learn policy $\psi^{(*)}$ for proposed parameter $\theta^{(*)}$ at
inverse temperature $\lambda_{i}$.

(4c) Run controlled SMC (Algorithm 1) with policy $\psi^{(*)}$ at proposed
parameter $\theta^{(*)}$ and inverse temperature $\lambda_{i}$ to obtain
likelihood estimator $\hat{p}(y_{1:T}|\theta^{(*)},\lambda_{i})$ and
trajectory $s_{0:T}^{(*)}$.

(4d) Update parameter $\theta^{(p)}=\theta^{(*)}$, trajectory
$s_{0:T}^{(p)}=s_{0:T}^{(*)}$ and likelihood estimator
$\hat{p}(y_{1:T}|\theta^{(p)},\lambda_{i})=\hat{p}(y_{1:T}|\theta^{(*)},\lambda_{i})$
with probability $\alpha(\theta^{(*)}|\theta^{(p)},\lambda_{i})$ in Equation
(34).

(5) If $\lambda_{i}=1$, set number of iterations as $I=i$ and terminate
iteration.
(6) Compute the model evidence estimator $\hat{p}(y_{1:T})$ as in Equation
(35).

Output: parameters $(\theta^{(p)})_{p=1}^{P}$, trajectories
$(s_{0:T}^{(p)})_{p=1}^{P}$ and model evidence estimator $\hat{p}(y_{1:T})$.

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

$\displaystyle\mathrm{ESS}(\Omega_{i}^{(1)}(\lambda),\ldots,\Omega_{i}^{(P)}(\lambda))=1{\Big{/}}\sum_{p=1}^{P}\Omega_{i}^{(p)}(\lambda)^{2},$
(31)

which is based on the unnormalized weights

$\displaystyle\omega_{i}^{(p)}(\lambda)=\omega(\theta^{(p)},s_{0:T}^{(p)}|\lambda_{i-1},\lambda)=\prod_{t=1}^{T}g_{\theta^{(p)}}(y_{t}|s_{t-1}^{(p)},s_{t}^{(p)})^{\lambda-\lambda_{i-1}}$
(32)

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:

$\displaystyle\lambda_{i}=\inf\left\\{\lambda\in(\lambda_{i-1},1]:\mathrm{ESS}\left(\Omega_{i}^{(1)}(\lambda),\ldots,\Omega_{i}^{(P)}(\lambda)\right)=\kappa_{\mathrm{ESS}}P\right\\},$
(33)

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 Dai et al.
(2020, Section 3) 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
Svensson, Schön, and Lindsten (2018), which requires storing the entire SMC
output for each parameter $\theta^{(p)}$, the unnormalized weight in Equation
(32) 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 (33) 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 5 in
Appendix B) 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 (Andrieu, Doucet, and Holenstein, 2010), which is necessary for the
validity of our approach. The use of conditional SMC within SMC2 was also
considered by Chopin, Jacob, and Papaspiliopoulos (2013), 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 Duan and Fulop (2015), 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 (29), 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 SMC2
approach. We refer readers to Dai et al. (2020, Section 2.2) 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 (Andrieu, Doucet, and Holenstein,
2010) 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 (Doucet et al., 2015; Sherlock et al., 2015), 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

$\alpha(\theta^{(*)}|\theta^{(p)},\lambda_{i})=\min\left\\{1,~{}\frac{p(\theta^{(*)})\hat{p}(y_{1:T}|\theta^{(*)},\lambda_{i})h_{i}(\theta^{(p)}|\theta^{(*)})}{p(\theta^{(p)})\hat{p}(y_{1:T}|\theta^{(p)},\lambda_{i})h_{i}(\theta^{(*)}|\theta^{(p)})}\right\\}.$
(34)

When the desired inverse temperature of $\lambda=1$ is reached, we terminate
the SMC2 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

$\hat{p}(y_{1:T})=\prod_{i=1}^{I}\left(\frac{1}{P}\sum_{p=1}^{P}\omega_{i}^{(p)}\right).$
(35)

Since estimators of the model evidence cannot be easily obtained with a single
PMMH chain targeting Equation (29), this is another strength of our SMC2
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$.

###  Consistency properties of adaptive SMC2 estimators

This section concerns asymptotic properties of our estimators of expectations
under the posterior distribution in Equation (29) and the model evidence
$p(y_{1:T})$. In the following, the notations $\stackrel{{\scriptstyle
p.}}{{\rightarrow}}$ and $\stackrel{{\scriptstyle d.}}{{\rightarrow}}$ denote
convergence in probability and distribution, respectively. We will establish
the weak law of large numbers (WLLN)

$\displaystyle\hat{\pi}(\varphi)\stackrel{{\scriptstyle
p.}}{{\longrightarrow}}\pi(\varphi),\quad\hat{p}(y_{1:T})\stackrel{{\scriptstyle
p.}}{{\longrightarrow}}p(y_{1:T}),$ (36)

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)

$\displaystyle\sqrt{P}(\hat{\pi}(\varphi)-\pi(\varphi))\stackrel{{\scriptstyle
d.}}{{\longrightarrow}}\mathcal{N}(0,\sigma^{2}(\varphi)),\quad\sqrt{P}(\hat{p}(y_{1:T})-p(y_{1:T}))\stackrel{{\scriptstyle
d.}}{{\longrightarrow}}\mathcal{N}(0,\sigma^{2}),$ (37)

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 SMC2 algorithm, it is clear how they can be
derived in our proofs, provided in the Appendix C.

We first consider Algorithm 4 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 SMC2
algorithm without the added complication of adaptation.

######  Theorem 2.

The estimators generated by SMC2 in Algorithm 4 without adaptation in Steps 2
and 4 satisfy the WLLNs and CLTs in Equations (36) and (37) for any bounded
and measurable function
$\varphi:\Theta\times\mathbb{S}^{T+1}\rightarrow\mathbb{R}$.

To study the adaptive SMC2 algorithm, we first note that the adaptation scheme
in Step 2 and Equation (33) should be seen as a finite sample approximation of
a limiting deterministic inverse temperature schedule $(\lambda_{i}^{*})$,
defined by the following scheme

$\displaystyle\lambda_{i}^{*}=\inf\left\\{\lambda\in(\lambda_{i-1}^{*},1]:\chi^{2}(p(d\theta,ds_{0:T}|y_{1:T},\lambda)~{}|~{}p(d\theta,ds_{0:T}|y_{1:T},\lambda_{i-1}^{*}))=\kappa_{\mathrm{ESS}}^{-1}-1\right\\},$
(38)

initialized at $\lambda_{0}^{*}=0$, where the above $\chi^{2}$-divergence

$\displaystyle\chi^{2}(p(d\theta,ds_{0:T}|y_{1:T},\lambda)~{}|~{}p(d\theta,ds_{0:T}|y_{1:T},\lambda_{i-1}^{*}))$
(39)
$\displaystyle=\frac{\int_{\Theta\times\mathbb{S}^{T+1}}\omega(\theta,s_{0:T}|\lambda_{i-1}^{*},\lambda)^{2}p(d\theta,ds_{0:T}|y_{1:T},\lambda_{i-1}^{*})}{(\int_{\Theta\times\mathbb{S}^{T+1}}\omega(\theta,s_{0:T}|\lambda_{i-1}^{*},\lambda)p(d\theta,ds_{0:T}|y_{1:T},\lambda_{i-1}^{*}))^{2}}-1$

depends on the importance weight $\omega$ defined in Equation (32). Next, we
formalize the adaptive nature of Step 4 using the framework of Beskos et al.
(2016). 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}^{P}S_{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 4, 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})$.

######  Assumption 3.

For time $t=1,\ldots,T$ and iteration $i\geq 1$, the observation density
$g_{\theta}(y_{t}|s_{t-1},s_{t})$, summary statistic $S_{i}$, importance
weight $\omega(\theta,s_{0:T}|\lambda_{i-1},\lambda_{i})$ and the non-negative
kernel $Q_{i}$ satisfy:

  1. (i)

$(\theta,s_{t-1},s_{t})\mapsto\log g_{\theta}(y_{t}|s_{t-1},s_{t})$ is bounded
on $\Theta\times\mathbb{S}\times\mathbb{S}$ for any $y_{t}\in\mathbb{Y}$;

  2. (ii)

$S_{i}:\Theta\rightarrow\mathbb{R}^{s}$ is bounded on $\Theta$;

  3. (iii)

$(\theta,s_{0:T},\lambda_{i-1},\lambda_{i})\mapsto\omega(\theta,s_{0:T}|\lambda_{i-1},\lambda_{i})$
is continuous at $(\lambda_{i-1}^{*},\lambda_{i}^{*})$ uniformly on
$\Theta\times\mathbb{S}^{T+1}$;

  4. (iv)

$(\theta,s_{0:T},\zeta_{i})\mapsto Q_{i}(\varphi|\theta,s_{0:T})$ is
continuous at $\zeta_{i}^{*}$ uniformly on $\Theta\times\mathbb{S}^{T+1}$ for
any bounded and measurable function
$\varphi:\Theta\times\mathbb{S}^{T+1}\rightarrow\mathbb{R}$;

  5. (v)

$(\theta,s_{0:T},\zeta_{i})\mapsto\nabla_{\zeta_{i}}Q_{i}(\varphi|\theta,s_{0:T})$
is well-defined, bounded and continuous at $\zeta_{i}^{*}$ uniformly on
$\Theta\times\mathbb{S}^{T+1}$ for any function
$\varphi:\Theta\times\mathbb{S}^{T+1}\rightarrow\mathbb{R}$.

######  Theorem 3.

Under Assumption 3, the estimators generated by adaptive SMC2 in Algorithm 4
satisfy the WLLNs and CLTs in Equations (36) and (37) for any bounded and
measurable function
$\varphi:\Theta\times\mathbb{S}^{T+1}\rightarrow\mathbb{R}$.

The crux of our arguments to establish Theorems 2 and 3 is to first cast the
rather involved SMC2 algorithm, without and with adaptation, as a particular
SMC sampler and adaptive SMC sampler in the frameworks of Del Moral et al.
(2006) and Beskos et al. (2016) that operate on an extended space. This
enables us to then invoke convergence results for standard SMC methods (Del
Moral, 2004; Chopin, 2004) and apply them to our SMC2 estimators. Such an
approach that exploits the specific properties of Algorithm 4 is not
applicable to the SMC2 method of Duan and Fulop (2015), which is based on the
tempered likelihood.

##  Applications

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 Woodford (2003), An and Schorfheide (2007), and Herbst and
Schorfheide (2016). This model has become a benchmark specification for the
analysis of monetary policy. Variants of this model have also been studied in
Aruoba, Cuba-Borda, and Schorfheide (2018) and Aruoba et al. (2021).

In the Appendix E, 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 (see, Bansal and Yaron, 2004; Bansal, Kiku, and Yaron,
2012a). This model has also been studied by Fulop et al. (2021).

###  Model setup

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

$\mathbb{E}_{t}\left[\sum_{s=0}^{\infty}\beta^{s}\left(\frac{(\mathsf{C}_{t+s}/\mathsf{A}_{t+s})^{1-\tau}-1}{1-\tau}+\chi_{M}\log\left(\frac{\mathsf{M}_{t+s}}{\mathsf{P}_{t+s}}\right)-\mathsf{H}_{t+s}\right)\right],$
(40)

subject to the budget constraint

$\mathsf{P}_{t}\mathsf{C}_{t}+\mathsf{B}_{t}+\mathsf{M}_{t}+\mathsf{T}_{t}=\mathsf{P}_{t}\mathsf{W}_{t}\mathsf{H}_{t}+\mathsf{R}_{t-1}\mathsf{B}_{t-1}+\mathsf{M}_{t-1}+\mathsf{P}_{t}\mathsf{D}_{t}+\mathsf{P}_{t}\mathsf{SC}_{t},$
(41)

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

$\mathsf{Y}_{t}(j)=\left(\frac{\mathsf{P}_{t}(j)}{\mathsf{P}_{t}}\right)^{-1/\nu}\mathsf{Y}_{t},$
(42)

and the price of the final good is

$\mathsf{P}_{t}=\left(\int_{0}^{1}\mathsf{P}_{t}(j)^{\frac{\nu-1}{\nu}}dj\right)^{\frac{\nu}{\nu-1}},$
(43)

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

$\mathsf{Y}_{t}(j)=\mathsf{A}_{t}\mathsf{N}_{t}(j),$ (44)

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

$\log\mathsf{A}_{t}=\log\gamma+\log\mathsf{A}_{t-1}+\log\mathsf{z}_{t},\quad\log\mathsf{z}_{t}=\rho_{z}\log\mathsf{z}_{t-1}+\varepsilon_{z,t}.$
(45)

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

$\mathsf{AC}_{t}(j)=\frac{\phi}{2}\left(\frac{\mathsf{P}_{t}(j)}{\mathsf{P}_{t-1}(j)}-\pi\right)^{2}\mathsf{Y}_{t}(j),$
(46)

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

$\mathbb{E}_{t}\left[\sum_{s=0}^{\infty}\beta^{s}\mathsf{Q}_{t+s|t}\left(\frac{\mathsf{P}_{t+s}(j)}{\mathsf{P}_{t+s}}\mathsf{Y}_{t+s}(j)-\mathsf{W}_{t+s}\mathsf{N}_{t+s}(j)-\mathsf{AC}_{t+s}(j)\right)\right],$
(47)

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

$\mathsf{G}_{t}=\left(1-\frac{1}{\mathsf{g}_{t}}\right)\mathsf{Y}_{t},$ (48)

where $\mathsf{g}_{t}$ is an exogenous process and is assumed to follow

$\log\mathsf{g}_{t}=(1-\rho_{g})\log\mathsf{g}+\rho_{g}\log\mathsf{g}_{t-1}+\varepsilon_{g,t},$
(49)

with $\varepsilon_{g,t}\sim\mathcal{N}(0,\sigma_{g}^{2})$.

A central bank sets the interest rate by an interest rate feedback rule

$\mathsf{R}_{t}=\mathsf{R}_{t}^{*1-\rho_{R}}\mathsf{R}_{t-1}^{\rho_{R}}\exp({\varepsilon_{R,t}}),$
(50)

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

$\mathsf{R}^{*}_{t}=\mathsf{r}\pi^{*}\left(\frac{\pi_{t}}{\pi^{*}}\right)^{\psi_{1}}\left(\frac{\mathsf{Y}_{t}}{\mathsf{Y}_{t}^{*}}\right)^{\psi_{2}},$
(51)

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

$\mathsf{P}_{t}\mathsf{G}_{t}+\mathsf{M}_{t-1}+\mathsf{R}_{t-1}\mathsf{B}_{t-1}=\mathsf{T}_{t}+\mathsf{M}_{t}+\mathsf{B}_{t}.$
(52)

_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

$\pi=\pi^{*},\quad\mathsf{r}=\gamma/\beta,\quad\mathsf{R}=\mathsf{r}\pi^{*},\quad\mathsf{c}=(1-\nu)^{1/\tau},\quad\mathsf{y}=\mathsf{g}\mathsf{c}.$
(53)

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

$\displaystyle\beta\mathbb{E}_{t}\Big{[}\exp(-\tau(\hat{\mathsf{c}}_{t+1}-\hat{\mathsf{c}}_{t})+\hat{\mathsf{R}}_{t}-\hat{\mathsf{z}}_{t+1}-\hat{\pi}_{t+1})\Big{]}=1,$
(54)
$\displaystyle\exp(\hat{\mathsf{c}}_{t}-\hat{\mathsf{y}}_{t})=\exp(-\hat{\mathsf{g}}_{t})-\frac{\phi\pi^{2}\mathsf{g}}{2}(\exp(\hat{\pi}_{t})-1)^{2},$
(55)
$\displaystyle\frac{1-\nu}{\nu\phi\pi^{2}}(\exp(\tau\hat{\mathsf{c}}_{t})-1)=(\exp(\hat{\pi}_{t})-1)\Big{[}(1-\frac{1}{2\nu})\exp(\hat{\pi}_{t})+\frac{1}{2\nu}\Big{]}$
$\displaystyle-\beta\mathbb{E}_{t}\Big{[}(\exp(\hat{\pi}_{t+1})-1)\exp(-\tau(\hat{\mathsf{c}}_{t+1}-\hat{\mathsf{c}}_{t})+\hat{\mathsf{y}}_{t+1}-\hat{\mathsf{y}}_{t}+\hat{\pi}_{t+1})\Big{]},$
(56)
$\displaystyle\hat{\mathsf{R}}_{t}=\rho_{R}\hat{\mathsf{R}}_{t-1}+(1-\rho_{R})\psi_{1}\hat{\pi}_{t}+(1-\rho_{R})\psi_{2}(\hat{\mathsf{y}}_{t}-\hat{\mathsf{g}}_{t})+\varepsilon_{R,t},$
(57)
$\displaystyle\hat{\mathsf{g}}_{t}=\rho_{g}\hat{\mathsf{g}}_{t-1}+\varepsilon_{g,t},$
(58)
$\displaystyle\hat{\mathsf{z}}_{t}=\rho_{z}\hat{\mathsf{z}}_{t-1}+\varepsilon_{z,t}.$
(59)

_Model solution._ The above equilibrium conditions form a non-linear rational
expectations system in the state variable

$s_{t}=(\hat{\mathsf{y}}_{t},\hat{\mathsf{c}}_{t},\hat{\pi}_{t},\hat{\mathsf{R}}_{t},\varepsilon_{R,t},\hat{\mathsf{g}}_{t},\hat{\mathsf{z}}_{t}),$
(60)

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 (1). In empirical studies, linear approximation methods (e.g., Sims,
2002) 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 (see, An and Schorfheide, 2007; Herbst and
Schorfheide, 2016).

However, as shown in Fernandez-Villaverde, Rubio-Ramirez, and Santos (2006),
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 (Schmitt-Grohe and Uribe, 2004, 2007;
Andreasen, Fernandez-Villaverde, and Rubio-Ramirez, 2018) 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 (1). Following An and
Schorfheide (2007) and Herbst and Schorfheide (2016), 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

$\displaystyle\textrm{YGR}_{t}$ $\displaystyle=$
$\displaystyle\gamma^{(Q)}+100(\hat{\mathsf{y}}_{t}-\hat{\mathsf{y}}_{t-1}+\hat{\mathsf{z}}_{t})+u_{1,t},$
(61) $\displaystyle\textrm{INF}_{t}$ $\displaystyle=$
$\displaystyle\pi^{(A)}+400\hat{\pi}_{t}+u_{2,t},$ (62)
$\displaystyle\textrm{INT}_{t}$ $\displaystyle=$
$\displaystyle\pi^{(A)}+\mathsf{r}^{(A)}+4\gamma^{(Q)}+400\hat{R}_{t}+u_{3,t},$
(63)

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

$\theta=(\tau,\nu,\kappa,1/\mathsf{g},\psi_{1},\psi_{2},\mathsf{r}^{(A)},\pi^{(A)},\gamma^{(Q)},\rho_{R},\rho_{g},\rho_{z},\sigma_{R},\sigma_{g},\sigma_{z},\sigma_{e,y},\sigma_{e,\pi},\sigma_{e,R}),$
(64)

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

$\beta=\frac{1}{1+\mathsf{r}^{(A)}/400},\quad\pi=1+\frac{\pi^{(A)}}{400},\quad\gamma=1+\frac{\gamma^{(Q)}}{100}.$
(65)

###  Implementation details

_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

$X_{\theta}(s)=c(\theta)+L(\theta)s+Q_{\theta}(s),\quad s\in\mathbb{S}.$ (66)

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

$z_{t}=\rho(\theta)z_{t-1}+\Sigma(\theta)\varepsilon_{t},\quad\varepsilon_{t}\sim\mathcal{N}(0_{d_{z}},I_{d_{z}}),\quad
t=1,\ldots,T,$ (67)

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

$s_{t}=A(\theta)s_{t-1}+B(\theta)\varepsilon_{t}+c_{\theta}(s_{t-1},\varepsilon_{t}),\quad\varepsilon_{t}\sim\mathcal{N}(0_{d_{z}},I_{d_{z}}),$
(68)

for $t=1,\ldots,T$, with the initial state

$s_{0}=B(\theta)\varepsilon_{0}+c_{\theta}((0_{d_{x}},\Sigma(\theta)\varepsilon_{0}),\varepsilon_{0}),\quad\varepsilon_{0}\sim\mathcal{N}(0_{d_{z}},I_{d_{z}}).$
(69)

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 D.1. In the framework of
Equation (2), this corresponds to an initial distribution and a Markov
transition kernel that are partially degenerate

$\displaystyle\mu_{\theta}(ds_{0})$
$\displaystyle=\delta_{\Phi_{\theta}^{(0)}(\varepsilon_{0})}(ds_{0})\mathcal{N}(\varepsilon_{0};0_{d_{z}},I_{d_{z}})d\varepsilon_{0},$
(70) $\displaystyle f_{\theta}(ds_{t}|s_{t-1})$
$\displaystyle=\delta_{\Phi_{\theta}(s_{t-1},\varepsilon_{t})}(ds_{t})\mathcal{N}(\varepsilon_{t};0_{d_{z}},I_{d_{z}})d\varepsilon_{t},$

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 (68) and (69). The log-linearized model can be easily obtained by
simply suppressing the $c_{\theta}$ terms in Equations (68) and (69).

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

$\displaystyle
g_{\theta}(y_{t}|s_{t})=\mathcal{N}(y_{t};d(\theta)+E(\theta)s_{t},F(\theta)),$
(71)

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 (61), (62), and
(63).

_Annealed controlled SMC._ The uncontrolled SMC within AC-SMC is taken as the
BPF; see Section 3.1 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 (70), 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

$\displaystyle\mathbb{F}_{0}$
$\displaystyle=\left\\{\psi_{t}(\varepsilon_{0})=\exp(-Q_{0}(\varepsilon_{0};\beta_{0})):|\beta_{0}|^{2}\leq\xi\right\\},$
(72) $\displaystyle\mathbb{F}_{t}$
$\displaystyle=\left\\{\psi_{t}(s_{t-1},\varepsilon_{t})=\exp(-Q(s_{t-1},\varepsilon_{t};\beta_{t})):\beta_{t}=\Lambda_{\theta}(\tilde{\beta}_{t}),|\tilde{\beta}_{t}|^{2}\leq\xi\right\\},$

for $t=1,\ldots,T$, where $Q_{0}(z;\beta)=z^{\top}Az+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

$Q(z,z^{\prime};\beta)={z}^{\prime\top}Az^{\prime}+{z}^{\prime\top}b+{z}^{\prime\top}Cz+z^{\top}Dz+z^{\top}e+f,$
(73)

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 (72) 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 (72) is well-specified when log-
linearization approximations are used to solve for the equilibrium conditions;
see the Appendix D.2 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 2) 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 D.3 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 D.4 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 (3.2) and
Algorithm 2. Lastly, as the form of the functions in Equation (72) are closed
under multiplication, policy refinements can be performed by updating the
coefficients (see, Appendix D.5).

###  Simulations

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 An and Schorfheide (2007, Table 3). 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 1 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.

Table 1: Log-likelihood Estimates: Linear Model

| ME: 5% |  | ME: 10%
---|---|---|---
|  | BPF |  |  |  | BPF |
| KF | $2^{12}$ | $2^{13}$ | $2^{14}$ | AC-SMC |  | KF | $2^{12}$ | $2^{13}$ | $2^{14}$ | AC-SMC
Mean | -2091.4 | -2441.1 | -2287.4 | -2194.6 | -2091.4 |  | -2162.2 | -2217.7 | -2193.2 | -2179.6 | -2162.2
Var |  | 3061.2 | 1491.2 | 416.9 | $2.3\times 10^{-4}$ |  |  | 150.9 | 132.9 | 45.9 | $1.0\times 10^{-5}$
Time |  | 1.26 | 2.58 | 4.93 | 2.59 |  |  | 1.3 | 2.65 | 5.04 | 2.59
| ME: 15% |  | ME: 20%
|  | BPF |  |  |  | BPF |
| KF | $2^{12}$ | $2^{13}$ | $2^{14}$ | AC-SMC |  | KF | $2^{12}$ | $2^{13}$ | $2^{14}$ | AC-SMC
Mean | -2265.7 | -2297.3 | -2284.9 | -2276.9 | -2265.7 |  | -2311.2 | -2323.0 | -2317.9 | -2314.4 | -2311.2
Var |  | 80.7 | 51.3 | 29.8 | $3.4\times 10^{-6}$ |  |  | 30.8 | 17.7 | 7.3 | $9.9\times 10^{-7}$
Time |  | 1.32 | 2.66 | 5.05 | 2.54 |  |  | 1.35 | 2.97 | 5.3 | 2.7

_Note:_ We generate a sequence of quarterly observations of YGR, INF and INT
from the linear model for 500 time periods. The true model parameter values
are the same as those used in An and Schorfheide (2007, Table 3). For each
measurement error (ME) setting, we implement 100 independent repetitions of
BPF and AC-SMC, based on the same sequence of simulated data, and compute the
sample means and variances of their log-likelihood estimates. The true log-
likelihood value is computed using a Kalman filter (KF). 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. The runtime is measured in
seconds.

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

Table 2: Log-Likelihood Estimates: Non-linear Model

| ME: 5% |  | ME: 10%
---|---|---|---
|  | BPF |  |  |  | BPF |
| KF | $2^{14}$ | $2^{15}$ | $2^{16}$ | AC-SMC |  | KF | $2^{14}$ | $2^{15}$ | $2^{16}$ | AC-SMC
Mean | -2204.4 | -2192.7 | -2147.3 | -2120.7 | -2089.4 |  | -2212.4 | -2164.5 | -2159.5 | -2156.2 | -2152.2
Var |  | 570.1 | 325.2 | 182.7 | $6.5\times 10^{-2}$ |  |  | 34.4 | 14.6 | 10.8 | $3.2\times 10^{-4}$
Time |  | 5.12 | 13.17 | 30.68 | 20.58 |  |  | 5.03 | 13.23 | 30.97 | 20.42
| ME: 15% |  | ME: 20%
|  | BPF |  |  |  | BPF |
| KF | $2^{14}$ | $2^{15}$ | $2^{16}$ | AC-SMC |  | KF | $2^{14}$ | $2^{15}$ | $2^{16}$ | AC-SMC
Mean | -2257.3 | -2222.1 | -2219.1 | -2218.2 | -2216.9 |  | -2359.7 | -2347.7 | -2346.6 | -2345.5 | -2344.3
Var |  | 13.8 | 6.7 | 2.7 | $3.9\times 10^{-4}$ |  |  | 9.2 | 4.4 | 1.9 | $5.0\times 10^{-4}$
Time |  | 5.04 | 13.21 | 31.16 | 20.78 |  |  | 5.06 | 13.22 | 31.25 | 20.59

_Note:_ We generate a sequence of quarterly observations of YGR, INF and INT
from the non-linear model for 500 time periods. The true model parameter
values are the same as those used in An and Schorfheide (2007, Table 3). For
each measurement error (ME) setting, we implement 100 independent repetitions
of BPF and AC-SMC, based on the same sequence of simulated data, and compute
the sample means and variances of their log-likelihood estimates. 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 runtime is
measured in seconds.

###  Real data application

We now apply our adaptive SMC2 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 (see, e.g., An
and Schorfheide, 2007; Herbst and Schorfheide, 2016), 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 (61), (62),
and (63): 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 1.

Figure 1: Time Series Data

_Note:_ Time series data include quarter-to-quarter per capita GDP growth
rates (YGR), annualized quarter-to-quarter inflation rates (INF), and
annualized nominal interest rates (INT), all of which are in percentages. 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 regions shaded in grey correspond to NBER recession periods.

Our adaptive SMC2 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 3 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 SMC2 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 2 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 (Doucet et al., 2015; Sherlock
et al., 2015).

Figure 2: Acceptance Rates

The middle column of Table 3 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.

Table 3: Prior Distributions and Posterior Estimates of Parameters

| Priors | Fixed ME (20%) | Free ME
---|---|---|---
Parameter | Support | Distribution | Mean | Std | Mean | Std
$\tau$ | $(0,\infty)$ | $\mathcal{TN}(2.00,0.50)$ | 2.07 | 0.42 | 1.73 | 0.43
$\nu$ | $(0,1)$ | $\mathcal{U}(0,1)$ | 0.43 | 0.27 | 0.68 | 0.20
$\kappa$ | $(0,\infty)$ | $\mathcal{TN}(0.20,0.20)$ | 0.76 | 0.15 | 0.06 | 0.04
$1/\mathsf{g}$ | $(0,1)$ | $\mathcal{U}(0,1)$ | 0.28 | 0.27 | 0.53 | 0.21
$\psi_{1}$ | $(0,\infty)$ | $\mathcal{TN}(1.50,0.25)$ | 1.73 | 0.18 | 1.80 | 0.23
$\psi_{2}$ | $(0,\infty)$ | $\mathcal{TN}(0.50,0.25)$ | 0.44 | 0.22 | 0.67 | 0.23
$\mathsf{r}^{(A)}$ | $(0,\infty)$ | $\mathcal{TN}(0.80,0.50)$ | 0.41 | 0.22 | 0.51 | 0.21
$\pi^{(A)}$ | $(0,\infty)$ | $\mathcal{TN}(4.00,2.00)$ | 4.81 | 0.61 | 3.65 | 0.30
$\gamma^{(Q)}$ | $\mathbb{R}$ | $\mathcal{N}(0.40,0.20)$ | 0.79 | 0.10 | 0.68 | 0.09
$\rho_{R}$ | $(0,1)$ | $\mathcal{U}(0,1)$ | 0.83 | 0.03 | 0.64 | 0.05
$\rho_{g}$ | $(0,1)$ | $\mathcal{U}(0,1)$ | 0.98 | 0.02 | 0.91 | 0.16
$\rho_{z}$ | $(0,1)$ | $\mathcal{U}(0,1)$ | 0.99 | 0.01 | 0.97 | 0.01
100$\sigma_{R}$ | $(0,\infty)$ | $\mathcal{TN}(0.30,4.00)$ | 0.18 | 0.02 | 0.03 | 0.01
100$\sigma_{g}$ | $(0,\infty)$ | $\mathcal{TN}(0.40,4.00)$ | 0.62 | 0.05 | 0.02 | 0.01
100$\sigma_{z}$ | $(0,\infty)$ | $\mathcal{TN}(0.40,4.00)$ | 0.08 | 0.02 | 0.01 | 0.00
$\sigma_{e,y}$ | $(0,\infty)$ | $\mathcal{TN}(0.12,1.00)$ | 0.12 | —— | 0.50 | 0.05
$\sigma_{e,\pi}$ | $(0,\infty)$ | $\mathcal{TN}(0.30,1.00)$ | 0.30 | —— | 1.37 | 0.10
$\sigma_{e,R}$ | $(0,\infty)$ | $\mathcal{TN}(0.49,1.00)$ | 0.49 | —— | 0.13 | 0.03
Model Evidence |  |  | -442.8 | -407.6

_Note:_ This table details the prior distributions and posterior estimates of
model parameters. The left column provides the exact distributional form, the
support, and the hyper-parameters of the prior distribution for each model
parameter. $\mathcal{N}$ stands for the normal distribution, $\mathcal{TN}$
the truncated normal distribution, and $\mathcal{U}$ the uniform distribution.
The model is estimated using our adaptive SMC2 algorithm with $P=1,024$
parameter particles and $N=1,024$ state particles within AC-SMC. The middle
and right columns report the posterior means and standard deviations of model
parameters when the standard deviations of the measurement errors fixed at 20%
of the sample standard deviations and treated as free parameters,
respectively.

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 3, 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 3, 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.

Figure 3:  Selected Smoothed Latent States

_Note:_ This figure illustrates (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 at 20% of their sample analogue and treating them as free
parameters, respectively.

Using the output of our adaptive SMC2 algorithm, we can examine the smoothing
distribution of the latent states, such as the exogenous states which are
usually used in policy analysis. Figure 3 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.

##  Concluding remarks

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 SMC2 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 SMC2 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. Bengtsson et al. (2008)) will call for the use of more
efficient approaches, such as the one developed here.

## Acknowledgements

This work was funded by CY Initiative of Excellence (grant ”Investissements
d’Avenir” ANR-16-IDEX-0008).

## References

  * An and Schorfheide (2007) An, S. and F. Schorfheide (2007).  Bayesian analysis of DSGE models.  Econometric reviews 26, 113–172.
  * Andreasen et al. (2018) Andreasen, M., J. Fernandez-Villaverde, and J. Rubio-Ramirez (2018).  The pruned state-space system for nonlinear DSGE models: Theory and empirical applications.  Review of Economic Studies 85, 1–49.
  * Andrieu et al. (2010) Andrieu, C., A. Doucet, and R. Holenstein (2010).  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, 269–342.
  * Aruoba et al. (2021) Aruoba, S., P. Cuba-Borda, K. Higa-Flores, S. Villalvazo, and F. Schorfheide (2021).  Piecewise-linear approximations and filtering for DSGE models with occasionally-binding constraints.  Review of Economic Dynamics.
  * Aruoba et al. (2018) Aruoba, S., P. Cuba-Borda, and F. Schorfheide (2018).  Macroeconomic dynamics near the ZLB: A tale of two countries.  Review of Economic Studies 85, 87–118.
  * Bansal et al. (2007) Bansal, R., R. Gallant, and G. Tauchen (2007).  Rational pessimism, rational exuberance, and asset pricing models.  Review of Economic Studies 74, 1005–1033.
  * Bansal et al. (2012a) Bansal, R., D. Kiku, and A. Yaron (2012a).  An empirical evaluation of the long-run risks model for asset prices.  Critical Finance Review 1, 183–221.
  * Bansal et al. (2012b) Bansal, R., D. Kiku, and A. Yaron (2012b).  An empirical evaluation of the long-run risks model for asset prices.  Critical Finance Review 1(1), 183–221.
  * Bansal et al. (2016a) Bansal, R., D. Kiku, and A. Yaron (2016a).  Risks for the long run: Estimation with time aggregation.  Journal of Monetary Economics 82, 52–69.
  * Bansal et al. (2016b) Bansal, R., D. Kiku, and A. Yaron (2016b).  Risks for the long run: Estimation with time aggregation.  Journal of Monetary Economics 82, 52–69.
  * Bansal and Yaron (2004) Bansal, R. and A. Yaron (2004).  Risks for the long run: A potential resolution of asset pricing puzzles.  The Journal of Finance 59(4), 1481–1509.
  * Beeler and Campbell (2012) Beeler, J. and J. Y. Campbell (2012).  The long-run risks model and aggregate asset prices: An empirical assessment.  Critical Finance Review 1(1), 141–182.
  * Bengtsson et al. (2008) Bengtsson, T., P. Bickel, and B. Li (2008).  Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems.  In Probability and statistics: Essays in honor of David A. Freedman, pp. 316–334. Institute of Mathematical Statistics.
  * Beskos et al. (2016) Beskos, A., A. Jasra, N. Kantas, and A. Thiery (2016).  On the convergence of adaptive sequential Monte Carlo methods.  Annals of Applied Probability 26(2), 1111–1146.
  * Borovička and Stachurski (2020) Borovička, J. and J. Stachurski (2020).  Necessary and sufficient conditions for existence and uniqueness of recursive utilities.  Journal of Finance 75, 1457–1493.
  * Briers et al. (2010) Briers, M., A. Doucet, and S. Maskell (2010).  Smoothing algorithms for state–space models.  Annals of the Institute of Statistical Mathematics 62, 61\.
  * Campbell and Shiller (1988) Campbell, J. Y. and R. J. Shiller (1988).  Stock prices, earnings, and expected dividends.  The Journal of Finance 43(3), 661–676.
  * Carpenter et al. (1999) Carpenter, J., P. Clifford, and P. Fearnhead (1999).  Improved particle filter for nonlinear problems.  IEE Proceedings-Radar, Sonar and Navigation 146, 2–7.
  * Chatterjee and Diaconis (2018) Chatterjee, S. and P. Diaconis (2018).  The sample size required in importance sampling.  The Annals of Applied Probability 28, 1099–1135.
  * Chopin (2004) Chopin, N. (2004).  Central limit theorem for sequential Monte Carlo methods and its application to Bayesian inference.  Annals of Statistics 32(6), 2385–2411.
  * Chopin et al. (2013) Chopin, N., P. E. Jacob, and O. Papaspiliopoulos (2013).  SMC2: an efficient algorithm for sequential analysis of state space models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 75, 397–426.
  * Chopin and Papaspiliopoulos (2020) Chopin, N. and O. Papaspiliopoulos (2020).  An introduction to sequential Monte Carlo.  Springer.
  * Christiano et al. (2005) Christiano, L., M. Eichenbaum, and C. Evans (2005).  Nominal rigidities and the dynamic effects of a shock to monetary policy.  Journal of Political Economy 113, 1–45.
  * Creal (2017) Creal, D. D. (2017).  A class of non-Gaussian state space models with exact likelihood inference.  Journal of Business & Economic Statistics 35(4), 585–597.
  * Dai et al. (2020) Dai, C., J. Heng, P. E. Jacob, and N. Whiteley (2020).  An invitation to sequential Monte Carlo samplers.  arXiv preprint arXiv:2007.11936.
  * DeJong et al. (2013) DeJong, D. N., R. Liesenfeld, G. V. Moura, J.-F. Richard, and H. Dharmarajan (2013).  Efficient likelihood evaluation of state-space representations.  Review of Economic Studies 80, 538–567.
  * Del Moral (2004) Del Moral, P. (2004).  Feynman-Kac formulae.  In Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications, pp. 47–93. Springer.
  * Del Moral et al. (2006) Del Moral, P., A. Doucet, and A. Jasra (2006).  Sequential Monte Carlo samplers.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68, 411–436.
  * Doucet et al. (2001) Doucet, A., N. De Freitas, and N. Gordon (2001).  An introduction to sequential Monte Carlo methods.  In Sequential Monte Carlo methods in practice, pp. 3–14. Springer.
  * Doucet et al. (2000) Doucet, A., S. Godsill, and C. Andrieu (2000).  On sequential Monte Carlo sampling methods for Bayesian filtering.  Statistics and computing 10, 197–208.
  * Doucet et al. (2015) Doucet, A., M. K. Pitt, G. Deligiannidis, and R. Kohn (2015).  Efficient implementation of Markov chain Monte Carlo when using an unbiased likelihood estimator.  Biometrika 102, 295–313.
  * Duan and Fulop (2015) Duan, J.-C. and A. Fulop (2015).  Density-tempered marginalized sequential Monte Carlo samplers.  Journal of Business & Economic Statistics 33, 192–202.
  * Epstein and Zin (1989) Epstein, L. G. and S. E. Zin (1989).  Substitution, risk aversion, and the temporal behavior of consumption and asset returns: A theoretical framework.  Econometrica 57(4), 937–969.
  * Fernandez-Villaverde et al. (2006) Fernandez-Villaverde, J., J. Rubio-Ramirez, and M. Santos (2006).  Convergence properties of the likelihood of computed dynamic models.  Econometrica 74, 891–910.
  * Fernandez-Villaverde et al. (2016) Fernandez-Villaverde, J., J. Rubio-Ramirez, and F. Schorfheide (2016).  Solution and estimation methods for dsge models.  In J. B. Taylor and H. Uhlig (Eds.), Handbook of Macroeconomics, Chapter 9, pp. 527–724. North Holland: Elsevier B.V.
  * Fernandez-Villaverde and Rubio-Ramirez (2007) Fernandez-Villaverde, J. and J. F. Rubio-Ramirez (2007).  Estimating macroeconomic models: A likelihood approach.  Review of Economic Studies 74, 1059–1087.
  * Finke et al. (2020) Finke, A., A. Doucet, and A. M. Johansen (2020).  Limit theorems for sequential MCMC methods.  Advances in Applied Probability 52, 377–403.
  * Fulop et al. (2021) Fulop, A., J. Heng, J. Li, and H. Liu (2021).  Bayesian estimation of long-run risk models using sequential Monte Carlo.  Journal of Econometrics Forthcoming.
  * Fulop and Li (2013) Fulop, A. and J. Li (2013).  Efficient learning via simulation: A marginalized resample-move approach.  Journal of Econometrics 176, 146–161.
  * Gallant et al. (2019) Gallant, A., M. Jahan-Parvar, and H. Liu (2019).  Does smooth ambiguity matter for asset pricing?  Review of Financial Studies 32, 3617–3666.
  * Godsill and Clapp (2001) Godsill, S. and T. Clapp (2001).  Improvement strategies for Monte Carlo particle filters.  In Sequential Monte Carlo methods in practice, pp. 139–158. Springer.
  * Gordon et al. (1993) Gordon, N. J., D. J. Salmond, and A. F. Smith (1993).  Novel approach to nonlinear/non-Gaussian Bayesian state estimation.  In IEE proceedings F (radar and signal processing), Volume 140, pp. 107–113. IET.
  * Gourieroux and Jasiak (2006) Gourieroux, C. and J. Jasiak (2006).  Autoregressive gamma processes.  Journal of Forecasting 25(2), 129–152.
  * Guarniero et al. (2017) Guarniero, P., A. M. Johansen, and A. Lee (2017).  The iterated auxiliary particle filter.  Journal of the American Statistical Association 112, 1636–1647.
  * Heng et al. (2020) Heng, J., A. N. Bishop, G. Deligiannidis, and A. Doucet (2020).  Controlled sequential Monte Carlo.  Annals of Statistics 48, 2904–2929.
  * Herbst and Schorfheide (2016) Herbst, E. and F. Schorfheide (2016).  Bayesian Estimation of DSGE Models.  Princeton University press.
  * Herbst and Schorfheide (2019) Herbst, E. and F. Schorfheide (2019).  Tempered particle filtering.  Journal of Econometrics 210, 26–44.
  * Jacob et al. (2015) Jacob, P. E., L. M. Murray, and S. Rubenthaler (2015).  Path storage in the particle filter.  Statistics and Computing 25, 487–496.
  * Judd (1992) Judd, K. L. (1992).  Projection methods for solving aggregate growth models.  Journal of Economic theory 58(2), 410–452.
  * Kong et al. (1994) Kong, A., J. S. Liu, and W. H. Wong (1994).  Sequential imputations and Bayesian missing data problems.  Journal of the American statistical association 89, 278–288.
  * Kydland and Prescott (1982) Kydland, F. and E. Prescott (1982).  Time to build and aggregate fluctuations.  Econometrica 50, 1345–1370.
  * Lucas (1978) Lucas, R. (1978).  Asset prices in an exchange economy.  Econometrica 46, 1429–1446.
  * Pitt and Shephard (1999) Pitt, M. K. and N. Shephard (1999).  Filtering via simulation: Auxiliary particle filters.  Journal of the American statistical association 94, 590–599.
  * Pohl et al. (2018) Pohl, W., K. Schmedders, and O. Wilms (2018).  Higher order effects in asset pricing models with long-run risks.  The Journal of Finance 73(3), 1061–1111.
  * Richard and Zhang (2007) Richard, J.-F. and W. Zhang (2007).  Efficient high-dimensional importance sampling.  Journal of Econometrics 141, 1385–1411.
  * Scharth and Kohn (2016) Scharth, M. and R. Kohn (2016).  Particle efficient importance sampling.  Journal of Econometrics 190, 133–147.
  * Schmitt-Grohe and Uribe (2004) Schmitt-Grohe, S. and M. Uribe (2004).  Solving dynamic general equilibrium models using a second-order approximation to the policy function.  Journal of Economic Dynamics and Control 28, 755–775.
  * Schmitt-Grohe and Uribe (2007) Schmitt-Grohe, S. and M. Uribe (2007).  Optimal simple and implementable monetary and fiscal rules.  Journal of Monetary Economics 54, 1702–1725.
  * Schorfheide et al. (2018) Schorfheide, F., D. Song, and A. Yaron (2018).  Identifying long-run risks: A Bayesian mixed-frequency approach.  Econometrica 86(2), 617–654.
  * Sherlock et al. (2015) Sherlock, C., A. H. Thiery, G. O. Roberts, and J. S. Rosenthal (2015).  On the efficiency of pseudo-marginal random walk Metropolis algorithms.  The Annals of Statistics 43, 238–275.
  * Sims (2002) Sims, C. A. (2002).  Solving linear rational expectations models.  Computational economics 20, 1–20.
  * Smets and Wouters (2003) Smets, F. and R. Wouters (2003).  An estimated dynamic stochastic general equilibrium model of the euro area.  Journal of European Economic Association 1, 1123–1175.
  * Svensson et al. (2018) Svensson, A., T. B. Schön, and F. Lindsten (2018).  Learning of state-space models with highly informative observations: A tempered sequential Monte Carlo solution.  Mechanical systems and signal processing 104, 915–928.
  * Weil (1989) Weil, P. (1989).  The equity premium puzzle and the risk-free rate puzzle.  Journal of Monetary Economics 24(3), 401–421.
  * Woodford (2003) Woodford, M. (2003).  Interest and Prices.  Princeton University press.

Appendix

##  Appendix A Proofs for annealed controlled sequential Monte Carlo

In this section, we detail the proofs of the results in Section 3.5 on
annealed controlled sequential Monte Carlo.

###### Proof of Proposition 1.

Using the relation $\psi^{*}=\psi\cdot\phi^{*}$, which holds as $\phi^{*}$ is
the optimal refinement of $\psi$ (Heng et al., 2020, Proposition 1), and the
form of the smoothing distribution $p(ds_{0:T}|y_{1:T},\theta,\lambda)$ in
Equation (13), we can decompose the log-density as

$\displaystyle\log\left(\frac{p(s_{0:T}|y_{1:T},\theta,\lambda)}{q^{\psi\cdot\phi}(s_{0:T})}\right)$
$\displaystyle=\log\left(\frac{\phi_{0}^{*}(s_{0})}{\phi_{0}(s_{0})}\right)+\log\left(\frac{q_{0}(\psi_{0}\cdot\phi_{0})}{q_{0}(\psi_{0}\cdot\phi_{0}^{*})}\right)$
$\displaystyle+\sum_{t=1}^{T}\log\left(\frac{\phi_{t}^{*}(s_{t-1},s_{t})}{\phi_{t}(s_{t-1},s_{t})}\right)+\sum_{t=1}^{T}\log\left(\frac{q_{t}(\psi_{t}\cdot\phi_{t}|s_{t-1})}{q_{t}(\psi_{t}\cdot\phi_{t}^{*}|s_{t-1})}\right).$
(A.1)

To upper bound the expectations and conditional expectations in (A), we apply
the log-sum inequality to obtain

$\displaystyle\log\left(\frac{q_{0}(\psi_{0}\cdot\phi_{0})}{q_{0}(\psi_{0}\cdot\phi_{0}^{*})}\right)\leq\int_{\mathbb{S}}\frac{q_{0}(ds_{0})\psi_{0}(s_{0})\phi_{0}(s_{0})}{q_{0}(\psi_{0}\cdot\phi_{0})}\log\left(\frac{\phi_{0}(s_{0})}{\phi_{0}^{*}(s_{0})}\right)=q_{0}^{\psi\cdot\phi}(\log(\phi_{0}/\phi_{0}^{*}))$
(A.2)

and

$\displaystyle\log\left(\frac{q_{t}(\psi_{t}\cdot\phi_{t}|s_{t-1})}{q_{t}(\psi_{t}\cdot\phi_{t}^{*}|s_{t-1})}\right)$
$\displaystyle\leq\int_{\mathbb{S}}\frac{q_{t}(ds_{t}|s_{t-1})\psi_{t}(s_{t-1},s_{t})\phi_{t}(s_{t-1},s_{t})}{q_{t}(\psi_{t}\cdot\phi_{t}|s_{t-1})}\log\left(\frac{\phi_{t}(s_{t-1},s_{t})}{\phi_{t}^{*}(s_{t-1},s_{t})}\right)$
$\displaystyle=q_{t}^{\psi\cdot\phi}(\log(\phi_{t}/\phi_{t}^{*})|s_{t-1}).$
(A.3)

Hence the KL divergence from $q^{\psi\cdot\phi}(ds_{0:T})$ to
$p(ds_{0:T}|y_{1:T},\theta,\lambda)$ can be upper bounded by

$\displaystyle\mathrm{KL}\left(p(ds_{0:T}|y_{1:T},\theta,\lambda)~{}|~{}q^{\psi\cdot\phi}(ds_{0:T})\right)=\int_{\mathbb{S}^{T+1}}\log\left(\frac{p(s_{0:T}|y_{1:T},\theta,\lambda)}{q^{\psi\cdot\phi}(s_{0:T})}\right)p(ds_{0:T}|y_{1:T},\theta,\lambda)$
$\displaystyle\leq
q_{0}^{*}(\log(\phi_{0}^{*}/\phi_{0}))+q_{0}^{\psi\cdot\phi}(\log(\phi_{0}/\phi_{0}^{*}))$
$\displaystyle+\sum_{t=1}^{T}(\mu_{t-1}^{*}\times
q_{t}^{*})(\log(\phi_{t}^{*}/\phi_{t}))+\sum_{t=1}^{T}(\mu_{t-1}^{*}\times
q_{t}^{\psi\cdot\phi})(\log(\phi_{t}/\phi_{t}^{*})).$ (A.4)

Equation (24) follows by recalling that $\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$. ∎

###### Proof of Theorem 1.

We first establish the backward recursion for the sequence
$(\varepsilon_{t}^{*})_{t=0}^{T}$. By Assumption 2$(i)$ and Jensen’s
inequality, the difference at the terminal time $T$ satisfies

$\displaystyle\varepsilon_{T}^{*}\leq
C_{T}\nu_{T}^{\psi}(\log(\phi_{T}^{*}/\phi_{T}))\leq
C_{T}\|\log\phi_{T}-\log\phi_{T}^{*}\|_{L^{2}(\nu_{T}^{\psi})}.$ (A.5)

Noting that $\phi_{T}=P_{T}^{\psi}w_{T}^{\psi}$, $\phi_{T}^{*}=w_{T}^{\psi}$
and $w_{T}^{\psi}\in\mathbb{G}_{T}^{\psi}$, we have $\varepsilon_{T}^{*}\leq
C_{T}e_{T}^{\psi}$ by applying Assumption 1$(ii)$. For time $t=0,\ldots,T-1$,
we consider the decomposition

$\displaystyle\varepsilon_{t}^{*}=\xi_{t}^{*}(\log(\phi_{t}^{*}/\tilde{\phi}_{t}))+\xi_{t}^{*}(\log(\tilde{\phi}_{t}/\phi_{t})),$
(A.6)

with $\tilde{\phi}_{t}=w_{t}^{\psi}q_{t+1}^{\psi}(\phi_{t+1})$. By the log-sum
inequality, we have

$\displaystyle\log(\phi_{t}^{*}/\tilde{\phi}_{t})=\log\left(\frac{q_{t+1}^{\psi}(\phi_{t+1}^{*})}{q_{t+1}^{\psi}(\phi_{t+1})}\right)=\log\left(\frac{q_{t+1}(\psi_{t+1}\phi_{t+1}^{*})}{q_{t+1}(\psi_{t+1}\phi_{t+1})}\right)\leq
q_{t+1}^{*}(\log(\phi_{t+1}^{*}/\phi_{t+1})).$ (A.7)

Hence using the definition of the distributions $(\xi_{t}^{*})_{t=0}^{T}$, the
first term in Equation (A.6) satisfies