The Block-Correlated Pseudo Marginal Sampler for State Space Models

Abstract

Pseudo Marginal Metropolis-Hastings (PMMH) is a general approach to Bayesian inference when the likelihood is intractable, but can be estimated unbiasedly. Our article develops an efficient PMMH method that scales up better to higher dimensional state vectors than previous approaches. The improvement is achieved by the following innovations. First, a novel block version of PMMH that works with multiple particle filters is proposed. Second, the trimmed mean of the unbiased likelihood estimates of the multiple particle filters is used. Third, the article develops an efficient auxiliary disturbance particle filter, which is necessary when the bootstrap disturbance filter is inefficient, but the state transition density cannot be expressed in closed form. Fourth, a novel sorting algorithm, which is as effective as previous approaches but significantly faster than them, is developed to preserve the correlation between the logs of the likelihood estimates at the current and proposed parameter values. The performance of the sampler is investigated empirically by applying it to nonlinear Dynamic Stochastic General Equilibrium models with relatively high state dimensions and with intractable state transition densities and to multivariate GARCH diffusion-driven volatility in the mean models. Although we only apply the method to state space models, the approach will be useful in a wide range of applications such as large panel data models and stochastic differential equation models with mixed effects.

Keywords:

• Block PMMH
• Correlated PMMH
• Dynamic stochastic general equilibrium (DSGE) model
• Trimmed mean likelihood estimate

1 Introduction

Pseudo marginal Metropolis-Hastings (PMMH) (Andrieu and Roberts 2009) is a Bayesian inference method for statistical models having an intractable likelihood, when the likelihood can be estimated unbiasedly. Our article develops a PMMH sampling method for efficiently estimating the parameters of complex state space models. The method scales better than current approaches in the number of observations and latent states and can handle state transition densities that cannot be expressed in closed form; for example, many dynamic stochastic general equilibrium (DSGE) models, which are a popular class of macroeconomic time series state space models, do not have closed form transition densities when second order expansions are used.

A key issue in efficiently estimating statistical models using a PMMH approach is that the variance of the log of the estimated likelihood grows with the number of observations and the dimension of the latent states (Deligiannidis, Doucet, and Pitt 2018). Pitt et al. (2012) show that to obtain a balance between computational time and the mixing of the Markov chain Monte Carlo (MCMC) chain, the number of particles used in the particle filter should be such that the variance of the log of the estimated likelihood is in the range 1–3, depending on the efficiency of the proposal for θ. Pitt et al. (2012) also show that the efficiency of PMMH schemes deteriorates exponentially as that variance increases; we further note that in many complex statistical applications, it is computationally very expensive ensuring that the variance of the log of the estimated likelihood is within the required range.

Deligiannidis, Doucet, and Pitt (2018) propose a more efficient PMMH scheme, called the correlated pseudo-marginal (CPM) method, which correlates the random numbers in the (single) particle filters used to construct the estimated likelihoods at the current and proposed values of the parameters. This dependence induces a positive correlation between the estimated likelihoods and reduces the variance of the difference in the logs of the estimated likelihoods which appear in the Metropolis-Hastings (MH) acceptance ratio. They show that the CPM scales up better with the number of observations compared to the standard pseudo marginal method of Andrieu, Doucet, and Holenstein (2010) when the state dimension is small.

Tran et al. (2016) propose an alternative correlated pseudo marginal approach to that of Deligiannidis, Doucet, and Pitt, calling it the block pseudo marginal (BPM) method; the BPM divides the random numbers used in constructing likelihood estimators into blocks and then updates the parameters jointly with one randomly chosen block of the random numbers in each MCMC iteration; this induces a positive correlation between the numerator and denominator of the MH acceptance ratio, similarly to the CPM. They show that for large samples the correlation of the logs of the estimated likelihoods at the current and proposed values is close to $1-1/S$, where S is the number of blocks. However, they do not apply the BPM method for estimating time series state space models using the particle filter.

Our article introduces a new PMMH sampler, referred to as the multiple PMMH algorithm (MPM), which extends the CPM method of Deligiannidis, Doucet, and Pitt (2018) and the BPM method of Tran et al. (2016). The MPM sampler is innovative and addresses the inherent estimation issues as follows: (a) The likelihood estimator is a trimmed mean of unbiased likelihood estimators from multiple independent particle filters (PFs); these PFs can be run in parallel. This approach reduces the variance of the likelihood estimator compared to using a single particle filter. The algorithm is exact when there is no trimming and approximate otherwise, but our empirical results suggest that the approximation is very close to being exact. (b) The unknown parameters, and only the random numbers used in one of the PFs, are updated jointly. This is a novel block version of PMMH that works with multiple particle filters; see Section 2.2 for details. (c) An auxiliary disturbance PF (ADPF) algorithm is proposed to estimate the likelihood efficiently for state space models such as DSGE models, where the bootstrap filter is very inefficient and the state transition density does not have a closed form so that an auxiliary PF cannot be constructed; see Section 2.3 for details. (d) A novel sorting algorithm, which is as effective as previous approaches but is significantly faster than them, is proposed to maintain the correlation between the logs of the likelihood estimates at the current and proposed values. The standard resampling step in the particle filter introduces discontinuities that break down the correlation between the logs of the likelihood estimates at the current and proposed values, even when the current and proposed parameters are close. Section 3.2 shows that the MPM sampler with the proposed sorting algorithm can maintain the correlation between the logs of the estimated likelihoods for relatively high dimensional state space models. Note that the proposed PMMH sampler works with any particle filter algorithm, including the tempered particle filter of Herbst and Schorfheide (2019), the standard bootstrap filter of Gordon, Salmond, and Smith (1993), the disturbance particle filter of Murray, Jones, and Parslow (2013), and the proposed auxiliary disturbance filter. When the number of state dimensions is greater than the number of disturbance dimensions, a disturbance particle filter is more efficient than a state-based particle filter. Additionally, the disturbance particle filter is useful for state space models with intractable state transition densities.

We illustrate the MPM sampler empirically and compare its performance to the CPM method, using a standard linear Gaussian state space model, a multivariate GARCH diffusion-driven volatility in the mean model and two nonlinear Dynamic Stochastic General Equilibrium (DSGE) models, using simulated and real data. Our work in estimating nonlinear DSGE models is also related to Fernández-Villaverde and Rubio-Ramírez (2007) and Hall, Pitt, and Kohn (2014) who use standard PMMH methods. We note that the MPM sampler will also be useful for other complex statistical models, such as panel data models and stochastic differential equation mixed effects models. The rest of the article is organized as follows. Section 2 introduces the state space model and discusses the MPM sampler. Section 3 presents results from both simulated and real data. Section 4 concludes with a discussion of our major results and findings. The article has an online supplement containing some further technical and empirical results. We have also made available the Matlab code to run all the empirical examples at https://github.com/davidgun1000/Block-Correlated-Pseudo-Marginal.

2 Bayesian Inference for State Space Models

Notation

We use the colon notation for collections of variables, that is, $a_t^{r:s}:=(a_t^r,\dots ,a_t^s)$ and for $t\le u,a_{t:u}^{r:s}:=(a_t^{r:s},\dots ,a_u^{r:s})$. Suppose $\left\{(Z_t,Y_t),t\ge 0\right\}$ is a stochastic process, with parameter θ, where the Y_t are the observations and the Z_t are the latent state vectors; random variables are denoted by capital letters and their realizations by lower case letters. We consider the state space model with $p(z_0|\theta )$ the density of $Z_0,$ $p(z_t|z_{t-1},\theta )$ the density of Z_t given $Z_{0:t-1}=z_{0:t-1}$ for $t\ge 1$, and $p(y_t|z_t,\theta )$ is the density of Y_t given $Z_{0:t}=z_{0:t}$, $Y_{1:t-1}=y_{1:t-1}$.

2.1 Bayesian Inference

The aim of Bayesian inference is to obtain the posterior distribution of the model parameters θ and the latent states $z_{0:T}$, given the observations $y_{1:T}$ and a prior distribution $p(\theta )$; that is, (1) $$p\left(\theta ,z_{0:T}|y_{1:T}\right)=p\left(y_{1:T}|\theta ,z_{0:T}\right)p\left(z_{0:T}|\theta \right)p\left(\theta \right)/p\left(y_{1:T}\right),$$(1) where $p(y_{1:T})=\int \int p(y_{1:T}|\theta ,z_{0:T})p(z_{0:T}|\theta )p(\theta )d{z}_{0:T}d\theta$ is the marginal likelihood. The likelihood, $p(y_{1:T}|\theta )=\int p(y_{1:T}|\theta ,z_{0:T})p(z_{0:T}|\theta )d{z}_{0:T},$ can be calculated exactly using the Kalman filter for linear Gaussian state space models (LGSS) and for linear DSGE models so that posterior samples can be obtained using an MCMC sampling scheme. However, the likelihood can be estimated unbiasedly, but not computed exactly, for nonlinear and non-Gaussian state space models and the nonlinear DSGE models described in Section 3.4.

The bootstrap particle filter (Gordon, Salmond, and Smith 1993) provides an unbiased estimator of the likelihood for a general state space model. Andrieu and Roberts (2009) and Andrieu, Doucet, and Holenstein (2010) show that it is possible to use this unbiased estimator of the likelihood to carry out exact Bayesian inference for the model parameters. They call this MCMC approach pseudo marginal Metropolis-Hastings (PMMH).

Algorithm 1

The Multiple PMMH (MPM) algorithm

• Set the initial values of ${\theta }^{(0)}$ arbitrarily.

• Sample $u_s\sim N(0,I)$ for $s=1,\dots ,S$, and run S (disturbance) particle filters (Algorithm 2 in Section 2.3) to estimate the likelihood ${\overline{\widehat{p}}}_N(y|\theta ,\tilde{u})$ as the trimmed mean of ${\widehat{p}}_N(y|\theta ,u_s),s=1,\dots ,S$; a 0% is the mean and a 50% trimmed mean is the median.

• For each MCMC iteration p, for $p=1,\dots ,P$,

  • Sample ${\theta }^{\prime }$ from the proposal density $q({\theta }^{\prime }|\theta )$.

  • Choose index s with probability $1/S$, sample ${\eta }_u\sim N(0,I)$, and set $u_s^{\prime }={\rho }_u{u}_s+\sqrt{1-{\rho }_u^2}{\eta }_u$.

  • Run S (disturbance) particle filters (Algorithm 2 in Section 2.3) to compute the estimate of likelihood ${\overline{\widehat{p}}}_N(y|{\theta }^{\prime },{\tilde{u}}^{\prime })$.

  • With the probability in (4), set $\overline{\widehat{p}}{(y_{1:T}|\theta ,\tilde{u})}^{(p)}=\overline{\widehat{p}}(y_{1:T}|{\theta }^{\prime },{\tilde{u}}^{\prime })$, ${\tilde{u}}^{(p)}={\tilde{u}}^{\prime }$, and ${\theta }^{(p)}={\theta }^{\prime }$; otherwise, set $\overline{\widehat{p}}{(y_{1:T}|\theta ,\tilde{u})}^{(p)}=\overline{\widehat{p}}{(y_{1:T}|\theta ,\tilde{u})}^{(p-1)}$, ${\tilde{u}}^{(p)}={\tilde{u}}^{(p-1)}$, and ${\theta }^{(p)}={\theta }^{(p-1)}$.

The nonlinear (second-order) DSGE models considered in Section 3.4 belong to the class of general nonlinear state space models whose state transition density is difficult to work with or cannot be expressed in closed form. In such cases, it is useful to express the model in terms of the density of its latent disturbance variables as it can be expressed in closed form. The posterior in (1) becomes (2) $$p\left(\theta ,{ϵ}_{1:T},z_0|y_{1:T}\right)\propto \prod _{t=1}^{T}p\left(y_t|{ϵ}_{1:t},\theta ,z_0\right)p\left({ϵ}_t\right)p\left(z_0|\theta \right)p\left(\theta \right),$$(2) which Murray, Jones, and Parslow (2013) call the disturbance state-space model. The standard state space model can be recovered from the disturbance state space model by using the deterministic function $F(z_{t-1},{ϵ}_t;\theta )\to z_t$. This gives us a state trajectory $z_{0:T}$ from any sample $({ϵ}_{1:T},z_0)$. In the disturbance state-space model the target becomes the posterior distribution over the parameters θ, the latent noise variables ${ϵ}_{1:T}$ and the initial state z₀, rather than $(\theta ,z_{0:T})$.

2.2 Multiple PMMH (MPM)

This section discusses algorithm 1, which is the proposed multiple PMMH (MPM) method that uses multiple particle filters to obtain the estimate of the likelihood. The novel version of the block-correlated PMMH that works with multiple particle filters to induce a high correlation between successive logs of the estimated likelihoods is also discussed. Section S1 of the online supplement provides a review of the standard pseudo marginal Metropolis-Hastings method of Andrieu and Roberts (2009). Our article focuses on MCMC based on the parameters, but it is straightforward to also use it to obtain posterior inference for the unobserved states; see, for example, Andrieu, Doucet, and Holenstein (2010).

Let U_s consists of all the random variables required to compute the unbiased likelihood estimate ${\widehat{p}}_N(y|\theta ,u_s)$ for the sth particle filter, with $p(u_s)$ the density of u_s . We now define the joint target density of θ and $\tilde{u}=(u_1,\dots ,u_S)$ as (3) $$p\left(\theta ,\tilde{u}|y_{1:T}\right)\propto {\overline{\widehat{p}}}_N\left(y|\theta ,\tilde{u}\right)p\left(\theta \right)\prod _{s=1}^{S}p\left(u_s\right),$$(3) where ${\overline{\widehat{p}}}_N(y|\theta ,\tilde{u})$ is the likelihood estimates obtained from S particle filters discussed below and $p(\theta )$ is the prior of θ. We update the parameters θ jointly with a randomly selected block u_s in each MCMC iteration, with $\mathrm{\text{Pr}}(S=s)=1/S$ for any $s=1,\dots ,S$. The selected block u_s is updated using $u_s^{\prime }={\rho }_u{u}_s+\sqrt{1-{\rho }_u^2}{\eta }_u$, where ρ_u is the nonnegative correlation between the random numbers u_s and $u_s^{\prime }$ and η_u is a standard normal vector of the same length as u_s . This is similar to component-wise MCMC whose convergence is well established in the literature; see, for example, Johnson, Jones, and Neath (2013). Using this scheme, the acceptance probability is (4) $$\begin{array}{c}\alpha \mathrm{}\left(\theta ,\tilde{u};{\theta }^{\prime },{\tilde{u}}^{\prime }\right)\\ =\mathrm{\text{min}}\mathrm{}\left(\mathrm{}1,\frac{\begin{array}{c}{\overline{\widehat{p}}}_N\left(y|{\theta }^{\prime },{\tilde{u}}^{\prime }=\left(u_1,\dots ,u_{s-1},u_s^{\prime },u_{s+1},\dots ,u_S\right)\right)\hfill \\ p\left({\theta }^{\prime }\right)q\left(\theta |{\theta }^{\prime }\right)\hfill \end{array}}{\begin{array}{c}{\overline{\widehat{p}}}_N\left(y|\theta ,\tilde{u}=\left(u_1,\dots ,u_{s-1},u_s,u_{s+1},\dots ,u_S\right)\right)\hfill \\ p\left(\theta \right)q\left({\theta }^{\prime }|\theta \right)\hfill \end{array}}\right).\end{array}$$(4)

Suppose that S particle filters are run in parallel. It is clear that the average of the S unbiased likelihood estimates from the particle filters, ${\overline{\widehat{p}}}_N(y|\theta ,\tilde{u})={\sum }_{s=1}^S{\widehat{p}}_N(y|\theta ,u_s)/S$, is also an unbiased estimate of the likelihood (Sherlock, Thiery, and Lee 2017). We now discuss a trimmed mean approach to obtain a likelihood estimate from multiple particle filters, which we then use in the MPM algorithm. For example, the 20% trimmed mean averages the middle 60% of the likelihood values. Although the trimmed mean does not give an unbiased estimate of the likelihood, we show in Section 3 that the posterior based on the trimmed means is a very good approximation to the exact posterior.

Algorithm 2

The Auxiliary Disturbance Particle Filter

Input: $u_{ϵ,1:T}^{1:N},u_{A,1:t-1}^{1:N}$, θ and N

Output: ${ϵ}_{1:T}^{1:N},z_{1:T}^{1:N},a_{1:T-1}^{1:N}$, and ${\overline{w}}_{1:T}^{1:N}$

For t = 1

• (1a) Generate ${ϵ}_1^i$ from $m({ϵ}_1^i|u_{ϵ,1}^i,\theta )$ and set $z_1^i=F(z_0^i,{ϵ}_1^i;\theta )$ for $i=1,\dots ,N$.

• (1b) Compute the unnormalized weights $w_1^i$, for $i=1,\dots ,N$.

• (1c) Compute the normalized weights, ${\overline{w}}_1^i=\frac{w_1^i}{{\sum }_{j=1}^N{w}_1^j}$ for $i=1,\dots ,N$.

For $t\ge 2$

• (2a) Sort the state particles $z_{t-1}^i$ or disturbance particles ${ϵ}_{t-1}^i$ using the multidimensional Euclidean sorting method in Section 2.4 and obtain the sorted index ζ_i for $i=1,\dots ,N$, and the sorted state particles, disturbance particles, and weights ${\tilde{z}}_{t-1}^i=z_{t-1}^{{\zeta }_i}$, ${\widetilde{ϵ}}_{t-1}^i={ϵ}_{t-1}^{{\zeta }_i}$ and ${\tilde{\overline{w}}}_t^i={\overline{w}}_{t-1}^{{\zeta }_i}$, for $i=1,\dots ,N$.

• (2b) Obtain the ancestor indices based on the sorted state particles or disturbance particles ${\tilde{a}}_{t-1}^{1:N}$ using the correlated multinomial resampling in Algorithm 3.

• (2c) Obtain the ancestor indices based on original order of the particles $a_{t-1}^i$, for $i=1,\dots ,N$.

• (2d) Generate ${ϵ}_t^i$ from $m({ϵ}_t^i|u_{ϵ,t}^i,\theta )$ and set $z_t^i=F(z_{t-1}^{a_{t-1}^i},{ϵ}_t^i;\theta )$, for $i=1,\dots ,N$.

• (2e) Compute the unnormalized weights $w_t^i$, for $i=1,\dots ,N$.

• (2f) Compute the normalized weights ${\overline{w}}_t^i=\frac{w_t^i}{{\sum }_{j=1}^N{w}_t^j}$ for $i=1,\dots ,N$.

2.3 The Auxiliary Disturbance Particle Filter

This section discusses the auxiliary disturbance particle filter (ADPF) that we use to estimate the likelihoods in the MPM sampler described in Section 2.2. It is particularly useful for state space models when the state dimension is substantially bigger than the disturbance dimension, and the state transition density is intractable. Suppose that $z_1=F(z_0,{ϵ}_1;\theta )$, where z₀ is the initial state vector with density $p(z_0|\theta )$, and $z_t=F(z_{t-1},{ϵ}_t;\theta )(t\ge 2)$, where ϵ_t is a vector of latent noise variables with density $p({ϵ}_t).$ Murray, Jones, and Parslow (2013) express the standard state-space model in terms of the latent noise variables ${ϵ}_{1:T}$, and call $${ϵ}_t\sim p\left({ϵ}_t\right),y_t|{ϵ}_{1:t}\sim p\left(y_t|{ϵ}_{1:t},z_0,\theta \right)=p\left(y_t|z_t,\theta \right),$$ $$t=1,\dots ,T,$$ the disturbance state-space model. We note that the conditional distribution of y_t depends on all the latent error variables, ${ϵ}_{1:t}.$

Let u be the value of the random variable used to obtain the unbiased estimate of the likelihood. It has two components $u_{ϵ,1:T}^{1:N}$ and $u_{A,1:T-1}^{1:N}$; $u_{ϵ,t}^i$ is the random variable used to generate the particles ${ϵ}_t^i$ given θ. We can write (5) $$\begin{array}{c}{ϵ}_1^i\sim m\left({ϵ}_1^i|u_{ϵ,1}^i,\theta \right),z_1^i=F\left(z_0^i,{ϵ}_1^i;\theta \right)\text{ and}\\ {ϵ}_t^i\sim m\left({ϵ}_t^i|u_{ϵ,t}^i,\theta \right),z_t^i=F\left(z_{t-1}^{a_{t-1}^i},{ϵ}_t^i;\theta \right),t\ge 2,\end{array}$$(5) where $m({ϵ}_t^i|u_{ϵ,t}^i,\theta )$ is the proposal density to generate ${ϵ}_t^i$, and z₀ is the initial state vector with density $p(z_0|\theta )$. Denote the distribution of $u_{ϵ,t}^i$ as ${\psi }_{ϵt}(·)$. For $t\ge 2$, let $u_{A,t-1}$ be the vector of random variables used to generate the ancestor indices $a_{t-1}^{1:N}$ using the resampling scheme $\mathcal{M}(a_{t-1}^{1:N}|{\overline{w}}_{t-1}^{1:N},z_{t-1}^{1:N})$ and define ${\psi }_{\mathit{\text{At}}-1}(·)$ as the distribution of $u_{A,t-1}$. Common choices for ${\psi }_{ϵt}(·)$ and ${\psi }_{\mathit{\text{At}}-1}(·)$ are iid $N(0,1)$ and iid $U(0,1)$ random variables, respectively. The argument $a_{t-1}^{1:N}$ means that $Z_{t-1}^{A_{t-1}^i}=z_{t-1}^{a_{t-1}^i}$ is the ancestor of $Z_t^i=z_t^i$.

Algorithm 2 takes the number of particles N, the parameters θ, the random variables used to generate the disturbance particles $u_{ϵ,1:T}^{1:N}$, and the random variables used in the resampling steps $u_{A,t-1}^{1:N}$ as the inputs; it outputs the set of state particles $z_{1:T}^{1:N}$, disturbance particles ${ϵ}_{1:T}^{1:N}$, ancestor indices $a_{1:T-1}^{1:N}$, and the normalized weights ${\overline{w}}_{1:T}^{1:N}$. At t = 1, the disturbance particles ${ϵ}_1^{1:N}$ are obtained as a function of the random numbers $u_{ϵ,1}^{1:N}$ using (5) in step (1a) and the state particles are obtained from $z_1^i=F(z_0^i,{ϵ}_1^i;\theta )$, for $i=1,\dots ,N$; the weights for all particles are then computed in steps (1b) and (1c).

Step (2a) sorts the state or disturbance particles from smallest to largest using the multidimensional Euclidean sorting procedure described in Section 2.4 to obtain the sorted disturbance particles, sorted state particles and weights. Algorithm 3 resamples the particles using the correlated multinomial resampling to obtain the sorted ancestor indices ${\tilde{a}}_{t-1}^{1:N}$. Step (2d) generates the disturbance particles ${ϵ}_t^{1:N}$ as a function of the random numbers $u_{ϵ,t}^{1:N}$ using (5) and the state particles are obtained from $z_t^i=F(z_{t-1}^{a_{t-1}^i},{ϵ}_t^i;\theta )$, for $i=1,\dots ,N$; we then compute the weights for all particles in steps (2e) and (2f). The disturbance particle filter provides an unbiased estimate of the likelihood ${\widehat{p}}_N(y|\theta ,u):={\prod }_{t=1}^T(N^{-1}{\sum }_{i=1}^N{w}_t^i),$ where $w_t^i=p(y_t|z_t^i,\theta )p({ϵ}_t^i)/m({ϵ}_t^i|u_{ϵ,t}^i,\theta )\text{ for }t=1,\dots ,T.$

We use the defensive mixture proposal density (Hesterberg 1995) (6) $$\begin{array}{c}m\left({ϵ}_t|u_{ϵ,t},\theta \right)=\pi p\left({ϵ}_t|\theta \right)+(1-\pi )q\left({ϵ}_t|\theta ,y_{1:T}\right),\hspace{1em}\text{ with }\\ 0<\pi \ll 1,\end{array}$$(6) to generate ϵ_t in the auxiliary disturbance particle filter; here $u_{ϵ,t}$ is the vector of random variables used to generate the particles ϵ_t given θ. If the observation density $p(y_t|z_t,\theta )$ is bounded, then (6) guarantees that the weights are bounded in the disturbance particle filter algorithm. In practice, we take $q({ϵ}_t|\theta ,y_{1:T})=N({ϵ}_t|{\widehat{\mu }}_t,{\widehat{\Sigma }}_t)$ and set $\pi =0.05$ in all the examples in Section 3; see Algorithm 4. If π = 1 in (6), then $m({ϵ}_t|\theta )=p({ϵ}_t|\theta )$ is the bootstrap disturbance particle filter. However, the empirical performance of this filter is usually poor because the resulting likelihood estimate is too variable.

Algorithm 3

Correlated Multinomial Resampling Algorithm

Input: $u_{A,t-1}$, sorted states ${\tilde{z}}_{t-1}^{1:N}$, sorted disturbances ${\widetilde{ϵ}}_{t-1}^{1:N}$, and sorted weights ${\tilde{\overline{w}}}_{t-1}^{1:N}$

Output: ${\tilde{a}}_{t-1}^{1:N}$

1. Compute the cumulative weights ${\widehat{F}}_{t-1}^N(j)={\sum }_{i=1}^j{\tilde{\overline{w}}}_{t-1}^i$ based on the sorted state particles $\left\{{\tilde{z}}_{t-1}^{1:N},{\tilde{\overline{w}}}_{t-1}^{1:N}\right\}$ or the sorted disturbance particles $\left\{{\widetilde{ϵ}}_{t-1}^{1:N},{\tilde{\overline{w}}}_{t-1}^{1:N}\right\}$

2. set ${\tilde{a}}_{t-1}^i=\underset{j}{\mathrm{\text{min}}}{\widehat{F}}_{t-1}^N(j)\ge u_{A,t-1}^i$ for $i=1,\dots ,N$. Note that ${\tilde{a}}_{t-1}^i$ for $i=1,\dots ,N$ is the ancestor index based on the sorted states or the sorted disturbances.

Algorithm 4

Constructing proposal for the auxiliary disturbance particle filter (ADPF)

Input: the number of particle filters S, the number of particles for each particle filter N, particles ${ϵ}_{1:S,1:T}^{1:N}$, weights ${\overline{w}}_{1:S,1:T}^{1:N}$, and ancestor indices $a_{1:S,1:T-1}^{1:N}$.

Output: the mean ${\widehat{\mu }}_t$ and the covariance matrix ${\widehat{\Sigma }}_t$, for $t=1,\dots ,T$.

1. Given the particles ${ϵ}_{1:S,1:T}^{1:N}$ with weights ${\overline{w}}_{1:S,1:T}^{1:N}$ and ancestor indices $a_{1:S,1:T-1}^{1:N}$ from the output of the disturbance particle filters, the ancestral tracing algorithm of Kitagawa (1996) is used to sample from particle approximations from the smoothing distribution $p({ϵ}_{1:T}|\theta ,y_{1:T})$. This consists of sampling one particle trajectory from each of the S particle filters in parallel. For each particle filter, we first sample an index j_s with the probability ${\overline{w}}_{s,T}^{j_s}$, tracing back its ancestral lineage $b_{s,1:T}^{j_s}(b_{s,T}^{j_s}=j_s\text{ and }{b}_{s,t-1}^{j_s}=a_{s,t-1}^{b_{s,t}^{j_s}})$ and choosing the particle trajectory ${ϵ}_{s,1:T}^{j_s}=({ϵ}_{s,1}^{b_{s,1}^{j_s}},\dots ,{ϵ}_{s,T}^{b_{s,T}^{j_s}})$ for $s=1,\dots ,S$.

2. The mean ${\widehat{\mu }}_t$ and the covariance matrix ${\widehat{\Sigma }}_t$ are

   ${\widehat{\mu }}_t:=\frac{1}{S}{\sum }_{s=1}^S{ϵ}_{s,t}^{b_{s,t}^{j_s}},\hspace{1em}\text{and}\hspace{1em}{\widehat{\Sigma }}_t:=\frac{1}{S-1}{\sum }_{s=1}^S({ϵ}_{s,t}^{b_{s,t}^{j_s}}-{\widehat{\mu }}_t){({ϵ}_{s,t}^{b_{s,t}^{j_s}}-{\widehat{\mu }}_t)}^{\top }.$

Algorithm 4 outlines how ${\widehat{\mu }}_t$ and the covariance matrix ${\widehat{\Sigma }}_t$ are obtained, for $t=1,\dots ,T$, at each iteration of the MPM algorithm. It is computationally cheap to obtain ${\widehat{\mu }}_t$ and ${\widehat{\Sigma }}_t$ for $t=1,\dots ,T$ because they are available from the output of the S disturbance particle filters run at each MCMC iteration and the ancestral tracing method is fast. The proposal mean ${\widehat{\mu }}_t$ and the covariance matrix ${\widehat{\Sigma }}_t$, for $t=1,\dots ,T$, obtained at iteration p are used to estimate the likelihood at the next iteration. At the first iteration, we use the bootstrap filter to initialize the mean ${\widehat{\mu }}_t$ and the covariance matrix ${\widehat{\Sigma }}_t$, for $t=1,\dots ,T$. Algorithm 5 gives the MPM algorithm with ADPF.

2.4 Multidimensional Euclidean Sorting Algorithm

It is important that the logs of the likelihood estimates evaluated at the current and proposed values of θ and $\tilde{u}$ are highly correlated to reduce the variance of $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{\overline{\widehat{p}}}_N(y|{\theta }^{\prime },{\tilde{u}}^{\prime })-\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{\overline{\widehat{p}}}_N(y|\theta ,\tilde{u})$ because this helps the Markov chain to mix well. However, the standard resampling step in the particle filter introduces discontinuities which breaks down the correlation between the logs of the likelihood estimates at the current and proposed values even when the current parameters θ and the proposed parameters ${\theta }^{\prime }$ are close. Sorting the particles from smallest to largest before resampling helps to preserve the correlation between the logs of the likelihood estimates at the current and proposed values (Deligiannidis, Doucet, and Pitt 2018). However, such simple sorting is unavailable for multidimensional state particles.

This section discusses Algorithm 6, which is the fast multidimensional Euclidean sorting algorithm used to sort the multidimensional state or disturbance particles in the particle filter algorithm described in Algorithm 2. We write the algorithm in terms of state particles, but similar steps can be implemented for the disturbance particles. The algorithm takes the particles $z_t^{1:N}$ and the normalized weights ${\overline{w}}_t^{1:N}$ as the inputs; it outputs the sorted particles ${\tilde{z}}_t^{1:N}$, sorted weights ${\tilde{\overline{w}}}_t^{1:N}$, and sorted indices ${\zeta }_{1:N}$. Let $z_t^i$ be the ith d-dimensional particles at time t, $z_t^i={(z_{t,1}^i,\dots ,z_{t,d}^i)}^{\top }$. Let $d(z_t^j,z_t^i)$ be the Euclidean distance between two multidimensional particles. The first step is to calculate the means of the multidimensional particles ${\overline{z}}_t^i=(1/d){\sum }_{k=1}^d{z}_{t,k}^i$ at time t for $i=1,\dots ,N$. The first sorting index, ζ₁, is the index of the multidimensional particle having the smallest value of ${\overline{z}}_t^i$ for $i=1,\dots ,N$. Therefore, the first sorted particle ${\tilde{z}}_t^1$ is $z_t^{{\zeta }_1}$ with its associated weight ${\tilde{\overline{w}}}_t^1={\overline{w}}_t^{{\zeta }_1}$. We then calculate the Euclidean distance between the selected multidimensional particle and the set of all remaining multidimensional particles, $d({\tilde{z}}_t^1,z_t^i)$, for all $i\ne {\zeta }_1$. The next step sorts the particles and weights according to the Euclidean distance from smallest to largest to obtain the sorted particles ${\tilde{z}}_t^{2:N}$, sorted weights ${\tilde{\overline{w}}}_t^{2:N}$, and sorted indices ${\zeta }_{2:N}$.

Algorithm 5

The Multiple PMMH (MPM) with ADPF algorithm

• Set the initial values of ${\theta }^{(0)}$ arbitrarily.

• Sample $u_s\sim N(0,I)$ for $s=1,\dots ,S$, run S (disturbance) particle filters to estimate the likelihood ${\overline{\widehat{p}}}_N(y|\theta ,\tilde{u})$ as the trimmed mean of ${\widehat{p}}_N(y|\theta ,u_s),s=1,\dots ,S$; a 0% trimmed mean is the mean and a 50% trimmed mean is the median, and run Algorithm 4 to construct the (initial) proposal for the auxiliary disturbance particle filter.

• For each MCMC iteration p, for $p=1,\dots ,P$,

  • Sample ${\theta }^{\prime }$ from the proposal density $q({\theta }^{\prime }|\theta )$.

  • Choose index s with probability $1/S$, sample ${\eta }_u\sim N(0,I)$, and set $u_s^{\prime }={\rho }_u{u}_s+\sqrt{1-{\rho }_u{}^2}{\eta }_u$.

  • Run S (disturbance) particle filters given in Algorithm 2 to compute the estimate of the likelihood ${\overline{\widehat{p}}}_N(y|{\theta }^{\prime },{\tilde{u}}^{\prime })$

  • Run Algorithm 4 to construct the proposal for the auxiliary disturbance particle filter at iteration p + 1.

  • With the probability in (4), set $\overline{\widehat{p}}{(y_{1:T}|\theta ,\tilde{u})}^{(p)}=\overline{\widehat{p}}(y_{1:T}|{\theta }^{\prime },{\tilde{u}}^{\prime }),{\tilde{u}}^{(p)}={\tilde{u}}^{\prime }$, and ${\theta }^{(p)}={\theta }^{\prime }$; otherwise, set $\overline{\widehat{p}}{(y_{1:T}|\theta ,\tilde{u})}^{(p)}=\overline{\widehat{p}}{(y_{1:T}|\theta ,\tilde{u})}^{(p-1)},{\tilde{u}}^{(p)}={\tilde{u}}^{(p-1)}$, and ${\theta }^{(p)}={\theta }^{(p-1)}$.

Algorithm 6 simplifies the sorting algorithm proposed by Choppala et al. (2016). In Choppala et al. (2016), the first sorting index is the index of the multidimensional particle having the smallest value along its first dimension. To select the second sorted index, ζ₂, we need to calculate the Euclidean distance between the first selected particle and the remaining multidimensional particles. Similarly, to select the third sorted index, ζ₃, we need to calculate the Euclidean distance between the second selected particle and the remaining multidimensional particles. This process is repeated N times, which is expensive for a large number of particles and a long time series. In contrast, Algorithm 6 only needs to calculate the Euclidean distance of the first selected particle and the remaining particles once. Deligiannidis, Doucet, and Pitt (2018) use the Hilbert sorting method of Skilling (2004) to order the multidimensional state particles, which requires calculating the Hilbert index for each multidimensional particle and is much more expensive than Algorithm 6.

Algorithm 6

Multidimensional Euclidean Sorting Algorithm

Input: $z_t^{1:N},{\overline{w}}_t^{1:N}$

Output: sorted particles ${\tilde{z}}_t^{1:N}$, sorted weights ${\tilde{\overline{w}}}_t^{1:N}$, sorted indices ${\zeta }_{1:N}$

• Calculate the means of the multidimensional particles ${\overline{z}}_t^i=\frac{1}{d}{\sum }_{k=1}^d{z}_{t,k}^i$ at time t for $i=1,\dots ,N$. The first sorting index, ζ₁, is the index of the multidimensional particle having the smallest value of ${\overline{z}}_t^i$ for $i=1,\dots ,N$. The first selected particle ${\tilde{z}}_t^1$ is $z_t^{{\zeta }_1}$ with its associated weight ${\tilde{\overline{w}}}_t^1={\overline{w}}_t^{{\zeta }_1}$.

• Calculate the Euclidean distance between the selected multidimensional particle and the set of all remaining multidimensional particles, $d({\tilde{z}}_t^1,z_t^i)$, for $\forall i\ne {\zeta }_1$.

• Sort the particles and weights according to the Euclidean distance from smallest to largest to obtain ζ_i for $i=2,\dots ,N$. The sorted particles ${\tilde{z}}_t^i=z_t^{{\zeta }_i}$ and ${\tilde{\overline{w}}}_t^i={\overline{w}}_t^{{\zeta }_i}$ for $i=2,\dots ,N$.

Table 1 shows the CPU time (in seconds) of the three sorting methods for different numbers of particles N and state dimensions d. From the table, the proposed sorting algorithm is much faster than that in Choppala et al. (2016) and the Hilbert sorting method. For example, for state dimensions d = 30 and N = 2000 particles, Algorithm 6 is 290 and 335 times faster than the Choppala et al. (2016) and the Hilbert sorting methods, respectively

Table 1 Comparing different sorting methods for various combinations of N and d in terms of CPU time (in seconds): I: The fast Euclidean sorting algorithm, II: Choppala et al. (2016) sorting method, and III: Hilbert sorting method.

d	N	I	II	III	I actual
10	500	1.00	30.50	87.50	0.0002
	1000	1.00	66.67	115.00	0.0003
	2000	1.00	172.40	141.00	0.0005
30	500	1.00	39.67	141.00	0.0003
	1000	1.00	113.00	214.00	0.0004
	2000	1.00	289.80	335.00	0.0005

NOTE: The entries in columns I–III are multiples of the time in column I. The column headed “I actual” is the actual time for column I; for example the 30.5 in column II and row 1 means that I is 30.5 times faster than the Choppala et al. (2016) sorting. The corresponding entry for I actual is 0.0002.

3 Examples

Section 3.1 discusses the inefficiency measures we use to compare the performance of different particle filters or PMMH samplers. Section 3.2 investigates empirically the performance of the proposed methods for estimating a high-dimensional linear Gaussian state space model. Section 3.3 discusses a multivariate GARCH diffusion-driven volatility in the mean model. Sections 3.4.1 and 3.4.2 apply the proposed samplers to estimate nonlinear small scale DSGE and medium scale DSGE models, respectively.

3.1 Definitions of Inefficiency

We use the inefficiency factor (IF) (also called the integrated autocorrelation time), $\text{IF}:=1+2{\sum }_{j=1}^{\infty }{\rho }_{\psi }(j)$, to measure the inefficiency of a PMMH sampler at estimating the posterior expectation of a univariate function $\psi (\theta )$ of θ; here, ${\rho }_{\psi }(j)$ is the jth autocorrelation of the iterates $\psi (\theta )$ in the MCMC chain after it has converged to its stationary distribution. We estimate the $\text{IF}$ using the CODA R package of Plummer et al. (2006). A low value of the $\text{IF}$ estimate suggests that the Markov chain mixes well. Our measure of the inefficiency of a PMMH sampler that takes computing time (CT) into account for a given parameter θ based on $\text{IF}$ is the time normalized inefficiency factor (TNIF) defined as $\text{TNIF}:=\text{IF}\times \text{CT}.$ For a given sampler, let ${\text{IF}}_{\mathrm{\text{MAX}}}$ and ${\text{IF}}_{\text{MEAN}}$ be the maximum and mean of the IF values over all the parameters in the model. The relative time normalized inefficiency factor (RTNIF) is a measure of the TNIF relative to the benchmark method, where the benchmark method depends on the example.

3.2 Linear Gaussian State Space Model

This section investigates (a) the efficiency of different approaches for obtaining likelihood estimates from S particle filters, and (b) the performance of different PMMH samplers for estimating a single parameter θ, regardless of the dimension of the states. We consider the linear Gaussian state space model discussed in Guarniero, Johansen, and Lee (2017) and Deligiannidis, Doucet, and Pitt (2018), where $\left\{X_t;t\ge 1\right\}$ and $\left\{Y_t;t\ge 1\right\}$ are ${\mathcal{R}}^d$ valued with $$\begin{array}{c}{Y}_t=X_t+W_t,\\ X_{t+1}=A_{\theta }{X}_t+V_{t+1},\end{array}$$ with $X_1\sim N(0_d,I_d),V_t\sim N(0_d,I_d)$, $W_t\sim N(0_d,I_d)$, and $A_{\theta }^{i,j}={\theta }^{|i-j|+1}$ for i < j; the true value of θ is 0.4. Although it is possible to use the more efficient fully adapted particle filter (Pitt et al. 2012), we use the bootstrap filter for all methods to show that the proposed methods are useful for models where it is difficult to use better particle filter algorithms.

The first study compares the 0%, 5%, 10%, 25%, and 50% trimmed means of the likelihood estimates obtained from S particle filters. The simulated data is generated from the model above with T = 200 and T = 300 time periods and d = 1, 5 and 10 dimensions.

Table 2 reports the sample variance of the log of the estimated likelihood obtained by using the five estimators of the likelihood for the d = 10 dimensional linear Gaussian state space model with T = 300 time periods. The table shows that: (a) there is no substantial reduction in the variance of the log of the estimated likelihood for the mean (0% trimmed mean). For example, the variance decreases by only a factor of 2.31 times when S increases from 20 to 1000 when N = 250. (b) The 5%, 10%, 25%, and 50% trimmed means of the individual likelihood estimates decrease the variance substantially as S and/or N increase. The 25% and 50% trimmed mean estimates have the smallest variance for all cases. Similar results are obtained for d = 5 with T = 200 and T = 300 and d = 10 with T = 200; see Tables S3, S4, and S5 in section S4 of the online supplement.

Table 2 Comparing the variance of the log of the estimated likelihood for five different estimators of the likelihood: I: Averaging the likelihood (0% trimmed mean), II: Averaging the likelihood (5% trimmed mean), III: Averaging the likelihood (10% trimmed mean), IV: Averaging the likelihood (25% trimmed mean), and V: Averaging the likelihood (50% trimmed mean), for d = 10 dimensional linear Gaussian state space model with T = 300.

N	S	I	II	III	IV	V	VI	V actual
100	1						563.65
	20	3.26	2.02	1.48	1.11	1.00		29.30
	100	11.47	2.60	1.80	1.22	1.00		5.73
	1000	70.45	2.32	1.58	1.01	1.00		0.66
250	1						629.83
	20	3.60	2.03	1.67	1.12	1.00		15.75
	100	11.90	2.48	1.56	1.12	1.00		3.20
	1000	76.55	2.24	1.62	1.12	1.00		0.32
1000	1						620.74
	20	3.69	1.94	1.51	0.99	1.00		6.14
	100	13.75	2.50	1.79	1.12	1.00		1.16
	1000	82.18	2.26	1.59	1.08	1.00		0.13

NOTE: The variance of the log of the estimated likelihood of a single particle filter is reported in column VI. The results are based on 1000 independent runs. The entries in columns headed I–V are relative to the variance in column V. The entries in column “V actual” are the actual variances for estimator V. For example, the entry on row 4, column I is 70.45, which means the variance of the log of the likelihood estimate is 70.45 times the corresponding variance in column V. The entries in column VI are relative to column V with S = 1000 particles.

Table 3 Comparing the performance of different PMMH samplers with different number of particle filters S and different number of particles N in each particle filter for estimating the linear Gaussian state space model using a simulate dataset with T = 300 and d = 10 dimensions.

	I	II		III		IV		V		V actual
N	50,000	250	500	250	500	250	500	250	500	250	500
S	1	100	100	100	100	100	100	100	100	100	100
$\widehat{\text{IF}}$	50.73	88.09	4.59	3.03	1.21	1.40	1.10	1.00	0.90	9.98	9.00
CT	8.18	1.00	2.00	1.00	2.00	1.00	2.00	1.00	2.00	0.99	1.98
$\widehat{\text{RTNIF}}$	415.08	88.10	9.17	3.03	2.42	1.40	2.20	1.00	1.80	9.88	17.82

NOTE: Sampler I: Correlated PMMH of Deligiannidis, Doucet, and Pitt (2018). Sampler II: MPM with 5% trimmed mean. Sampler III: MPM with 10% trimmed mean. Sampler IV: MPM with 25% trimmed mean. Sampler V: MPM with 50% trimmed mean. Time denotes the time taken in seconds for one iteration of the method. The $\widehat{\text{IF}}$, CT, and $\widehat{\text{RTNIF}}$ entries in columns headed I to V are relative to the entries in column V (MPM with 50% trimmed mean with N = 250 and S = 100). The entries in column “V actual” are the actual $\widehat{\text{IF}}$, CT, and $\widehat{\text{RTNIF}}$ values for estimator V. For example, the entry on row 5, column I is 415.08, which means the $\widehat{\text{RTNIF}}$ is 415.08 times the corresponding $\widehat{\text{RTNIF}}$ in column V (MPM with 50% trimmed mean with N = 250 and S = 100).

Table 4 Comparing the performance of different PMMH samplers with different number of particle filters S and different number of particles N in each particle filter for estimating multivariate GARCH diffusion-driven volatility in mean model using a real dataset with T = 100 and d = 30.

Param.	I			II
N	250	250	500	250	250
S	100	200	100	100	200
ψ	183.559	124.258	227.685	151.745	80.073
μ	74.070	49.250	384.767	74.485	45.554
a	46.159	45.662	126.968	32.721	14.753
${\tau }^2$	102.234	70.255	218.866	78.171	42.075
$\widehat{{\text{IF}}_{\mathrm{\text{MAX}}}}$	183.559	124.258	384.767	151.745	80.073
${\widehat{\text{TNIF}}}_{\mathrm{\text{MAX}}}$	523.143	877.262	3351.321	432.473	565.315
${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$	1	1.677	6.406	0.827	1.081
$\widehat{{\text{IF}}_{\text{MEAN}}}$	101.505	72.356	239.571	84.280	45.614
${\widehat{\text{TNIF}}}_{\text{MEAN}}$	289.289	510.833	2077.953	240.198	322.035
${\widehat{\text{RTNIF}}}_{\text{MEAN}}$	1	1.766	7.183	0.830	1.113
Time	2.85	7.06	8.71	2.85	7.06

NOTE: The model parameters are ψ, μ, a, and ${\tau }^2$. Sampler I: MPM with 25% trimmed mean. Sampler II: MPM with 50% trimmed mean. Time is the time in seconds for one iteration of the method. The RTNIF is the TNIF relative to the MPM with 25% trimmed mean with S = 100 and N = 250.

Table 5 Comparing the performance of different PMMH samplers using different number of particle filters S and different number of particles N in each particle filter for estimating the small scale DSGE model using a real dataset with T = 124 observations.

	I	II	III	IV	V	VI	VII	VIII	IX	VI actual
S	100	100	100	100	100	100	100	100	1	100
N	250	100	100	100	100	100	100	250	20000	100
${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}}$	3.73	0.25	0.22	0.25	0.28	1.00	0.39	NA	2.42	604.72
${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$	36.27	$\mathbf{1.26}$	$\mathbf{1.11}$	$\mathbf{1.25}$	1.54	$\mathbf{1.00}$	1.96	NA	114.37	66.52
${\widehat{\text{IF}}}_{\text{MEAN}}$	5.88	0.49	0.50	0.49	0.61	1.00	0.59	NA	3.32	190.74
${\widehat{\text{TNIF}}}_{\text{MEAN}}$	57.23	$\mathbf{2.47}$	$\mathbf{2.49}$	$\mathbf{2.44}$	3.41	$\mathbf{1.00}$	2.93	NA	157.02	20.98
Time	9.73	5.00	5.00	5.00	5.55	1.00	5.00	6.82	47.27	0.11

NOTE: Sampler I: MPM (0% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, ADPF). Sampler II: MPM (10% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, ADPF). Sampler III: MPM (25% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, ADPF). Sampler IV: MPM (25% trimmed mean, ${\rho }_u=0$, disturbance sorting, ADPF). Sampler V: MPM (25% trimmed mean, ${\rho }_u=0.9$, state sorting, ADPF). Sampler VI: DA-MPM (25% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, ADPF). Sampler VII: MPM (50% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, ADPF). Sampler VIII: MPM (0% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, bootstrap). Sampler IX: Correlated PMMH of Deligiannidis, Doucet, and Pitt (2018). Time is the time taken in seconds for one iteration of the method. The ${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{IF}}}_{\text{MEAN}}$, and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ entries in columns headed I to IX are relative to the entries in column VI. The entries in column “VI actual” are the actual ${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{IF}}}_{\text{MEAN}}$, and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ for estimator VI. The numbers in bold are for the sampler with similar values of the ${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$ and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$.

We now examine empirically the ability of the proposed MPM samplers to maintain the correlation between successive values of the log of the estimated likelihood for the 10 dimensional linear Gaussian state space model with T = 300. We ran the different MPM approaches for 1000 iterations holding the current parameter at $\theta =0.4$ and the proposed parameter at ${\theta }^{\prime }=0.4,0.399,0.385$. At each iteration we generated $\tilde{u}$, where $\tilde{u}=(u_1,\dots ,u_S)$, and ${\tilde{u}}^{\prime }$ and obtained $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{\overline{\widehat{p}}}_N(y|\theta ,{\tilde{u}}^{(s)})$ and $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{\overline{\widehat{p}}}_N(y|{\theta }^{\prime },{\tilde{u}}^{(s)\prime })$ for the MPM approaches and computed their sample correlations.

Figure S5 in section S4 of the online supplement reports the correlation estimates of the log of estimated likelihood obtained using different MPM approaches. The figure shows that: (a) when the current and proposed values of the parameters are equal to 0.4, all MPM methods maintain a high correlation between the logs of the estimated likelihoods. The estimated correlations between logs of the estimated likelihoods when the number of particle filters S is 100 are about 0.99; (b) when the current parameter θ and the proposed parameter ${\theta }^{\prime }$ are (slightly) different, the MPM methods can still maintain some of the correlations between the log of the estimated likelihoods. The MPM methods with 25% and 50% trimmed means of the likelihood perform the best to maintain the correlation between the logs of estimated likelihoods.

We now compare the efficiency of the different sampling schemes for estimating the parameter θ of the linear Gaussian state space model. In all the examples, we run the samplers for $25,000$ iterations, with the initial 5000 iterations discarded as warm up. The particle filter and the parameter samplers are implemented in Matlab running on a high performance computer cluster with 20 CPU cores. The optimal number of CPU cores required for the MPM algorithm is equal to the number of particle filters S, which means the properties of the sampler can be easily tuned to provide maximum parallel efficiency on a large range of hardware. The adaptive random walk proposal of Roberts and Rosenthal (2009) is used for all samplers.

Figure S2 in section S4 of the online supplement provides the trace plots of the parameter θ estimated using the MPM algorithm with the different trimmed means for estimating the 10 dimensional linear Gaussian state space model with T = 200. The figure shows that: (a) the MCMC chain gets stuck for the 5% trimmed mean approach with S = 100 and N = 250 and N = 500 and for the 10% trimmed mean approach with S = 100 and N = 100; (b) the 10% trimmed mean approach requires at least S = 100 and N = 250 to make the MCMC chain mix well; (c) the 25% trimmed mean approach is the most efficient sampler as its MCMC chain does not get stuck even with S = 100 and N = 100. Figure S1 in section S4 of the online supplement compares the MPM algorithms with and without blocking with S = 100 and N = 250 for estimating the 10 dimensional linear Gaussian state space model with T = 300. The figure shows that the MPM (no blocking) with 5% and 10% trimmed means get stuck and the MPM (no blocking) with 25% trimmed means mixes poorly. The $\widehat{\text{IF}}$ of the MPM (blocking) with 25% trimmed mean is 9.624 times smaller than the MPM (no blocking) with 25% trimmed mean. All three panels show that the MPM (blocking) algorithms using trimmed means perform better than the MPM (no blocking), demonstrating the usefulness of the MPM (blocking) sampling scheme.

Table 3 reports the $\widehat{\text{IF}}$, $\widehat{\text{TNIF}}$, and $\widehat{\text{RTNIF}}$ values for the parameter θ in the linear Gaussian state space model with d = 10 dimensions and T = 300 time periods estimated using the following five MCMC samplers: (a) the correlated PMMH of Deligiannidis, Doucet, and Pitt (2018), (b) the MPM with 5% trimmed mean, (c) the MPM with 10% trimmed mean, (d) the MPM with 25% trimmed mean, and (e) the MPM with 50% trimmed mean. In this article, instead of using the Hilbert sorting method, we implement the correlated PMMH using the fast multidimensional Euclidean sorting method in Section 2.4. The computing time reported in the table is the time to run a single particle filter for the CPM and S particle filters for the MPM approach. The table shows the following. (a) The correlated PMMH requires more than $50,000$ particles to improve the mixing of the MCMC chain for the parameter θ. (b) The CPU time for running a single particle filter with $N=50,000$ particles is 4.09 times higher than running multiple particle filters with S = 100 and N = 500. The MPM method can be much faster than the CPM method if it is run using high-performance computing with a large number of cores. (c) The MPM allows us to use a much smaller number of particles for each independent PF and these multiple PFs can be run independently. (d) The $\widehat{\text{IF}}$ values for the parameter θ estimated using the correlated PMMH with $N=50,000$ particles are approximately 11, 42, 46, and 56 times larger than the 5%, 10%, 25%, and 50% trimmed means approaches with S = 100 and N = 500 particles, respectively. (e) In terms of $\widehat{\text{RTNIF}}$, the 5%, 10%, 25%, and 50% trimmed means with S = 100 and N = 500 are approximately 45, 172, 189, and 231 times smaller than the correlated PMMH with $N=50,000$ particles. (6) The best sampler for this example is the 50% trimmed mean approach with S = 100 and N = 250 particles. Table S6 in section S4 of the online supplement gives similar results for the 10 dimensional linear Gaussian state space model with T = 200. Figure S4 in section S4 of the online supplement plots the kernel density estimates of the parameter θ estimated using Metropolis-Hastings algorithm with the (exact) Kalman filter method and the MPM algorithm with 5%, 10%, 25%, and 50% trimmed means of the estimated likelihood. The figure shows that the approximate posteriors obtained by various approaches using trimmed means are very close to the true posterior. Similar results are obtained for the case d = 10 dimensions and T = 200 time periods given in Figure S3 in section S4 of the online supplement. Figure S3 also shows that the approximate posteriors obtained using the MPM with the 5%, 10%, 25%, and 50% trimmed means of the likelihood are very close to the exact posterior obtained using the correlated PMMH.

In summary, the example suggests that: (a) for a high dimensional state space model, there is no substantial reduction in the variance of the log of the estimated likelihood for the method which uses the average of the likelihood from S particle filters when S and/or N increases. Methods that use trimmed means reduce the variance substantially when S and/or N increases; (b) the 25% and 50% trimmed mean approaches give the smallest variance of the log of the estimated likelihood for all cases; (c) the MPM method with 50% trimmed mean is best at maintaining the correlation between the logs of the estimated likelihoods in successive iterates; (d) the MPM (blocking) method is more efficient than the MPM (no blocking) for the same values of S and N; (e) all approximate approaches that use the trimmed means of the likelihood give accurate approximations to the true posterior; (f) the best sampler for estimating the 10 dimensional linear Gaussian state space model with T = 300 is the 50% trimmed mean approach with S = 100 and N = 250 particles; (g) the MPM allows us to use much smaller number of particles for each independent PF and these multiple PFs can be run independently. These approaches can be made much faster if they are run using a high-performance computer cluster with a large number of cores.

3.3 Multivariate GARCH Diffusion-Driven Volatility in Mean Model

This section investigates the performance of the proposed MPM samplers for estimating large multivariate GARCH diffusion-driven volatility in mean models, where each of the log volatility processes follows a GARCH (Generalized Autoregressive Conditional Heteroscedasticity) diffusion process (Shephard, Chib, and Pitt 2004). The GARCH diffusion model does not have a closed form state transition density, making it challenging to estimate using the Gibbs type MCMC sampler used in Cross et al. (2021). It is well-known that the Gibbs sampler is inefficient for generating parameters for a diffusion process, especially for the variance parameter ${\tau }^2$ (Stramer and Bognar 2011). Conversely, MPM samplers can efficiently estimate models with non-closed form state transition densities, making them well-suited for the GARCH diffusion model. This section demonstrates that MPM samplers can estimate this model efficiently and overcome the limitations of the Gibbs sampler.

Suppose that P_t is a $d\times 1$ vector of daily stock prices and define $y_t:=\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{P}_t-\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}{P}_{t-1}$ as the log-return of the stocks. Let $h_{i,t}$ be the log-volatility process of the ith stock at time t. We also define, $h_{\cdot ,t}:={(h_{1,t},\dots ,h_{d,t})}^{\top }$ and $h_{i,\cdot }:={(h_{i,1},\dots ,h_{i,T})}^{\top }$. The model for y_t is (7) $$y_t=A{h}_{·,t}+{ϵ}_t,{ϵ}_t\sim N\left(0,{\Sigma }_t\right),$$(7) where A is a d × d matrix that captures the effects of log-volatilities on the levels of the variables. The time-varying error covariance matrix Σ_t depends on the unobserved latent variables $h_{\cdot ,t}$ such that (8) $${\Sigma }_t:=\text{diag}\left(\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(h_{1,t}\right),\dots ,\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(h_{d,t}\right)\right).$$(8)

Each log-volatility $h_{i,t}$ is assumed to follow a continuous time GARCH diffusion process (Shephard, Chib, and Pitt 2004) satisfying (9) $$\text{ d}{h}_{i,t}=\left\{\alpha \left(\mu -\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(h_{i,t}\right)\right)\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(-h_{i,t}\right)-\frac{{\tau }^2}{2}\right\}\text{ d}t+\tau \text{ d}{W}_{i,t},$$(9) for $i=1,\dots d$, where the $W_{i,t}$ are independent Wiener processes. The following Euler scheme approximates the evolution of the log-volatilities $h_{i,t}$ in (9) by placing M–1 evenly spaced points between times t and t + 1. We denote the intermediate volatility components by $h_{i,t,1},\dots ,h_{i,t,M-1}$, and it is convenient to set $h_{i,t,0}=h_{i,t}$ and $h_{i,t,M}=h_{i,t+1}$. The equation for the Euler evolution, starting at $h_{i,t,0}$ is (see, e.g., Stramer and Bognar 2011, p. 234) is (10) $$\begin{array}{c}{h}_{i,t,j+1}|h_{i,t,j}\\ \sim N\left(h_{i,t,j}+\left\{\alpha \left(\mu -\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(h_{i,t,j}\right)\right)\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(-h_{i,t,j}\right)-\frac{{\tau }^2}{2}\right\}\delta ,{\tau }^2\delta \right),\end{array}$$(10) for $j=0,\dots ,M-1$, where $\delta =1/M$. The initial state of $h_{i,t}$ is assumed normally distributed N(0, 1) for $i=1,\dots ,d$. For this example, we assume for simplicity that the parameters μ, α, and ${\tau }^2$ are the same across all stocks and $A^{i,j}={\psi }^{|i-j|+1}$ for i < j. This gives the same number of parameters regardless of the dimension of the stock returns so that we can focus on comparing different PMMH samplers.

We first report on a simulation study for the above model. We simulated ten independent datasets from the model described above with d = 20 dimensions and T = 100 observations. The true parameters are α = 2, $\mu =1.81,{\tau }^2=0.38$, and $\psi =0.01$. The priors are $\alpha \sim G(1,1),{\tau }^2\sim G(0.5,0.5),\mu \sim G(1,1)$, and $\psi \sim N(0,1)$. These prior densities are non-informative and cover almost all possible values in practice. We use M = 3 latent points for the Euler approximations of the state transition densities. For each dataset, we run the MPM sampler with the 25% trimmed mean approach for $25,000$ iterations; the initial 5000 iterations are discarded as warm up. The particle filter and the parameter samplers are implemented in Matlab running on a high performance computer cluster with 20 cores. The adaptive random walk proposal of Roberts and Rosenthal (2009) is used for all samplers.

Figures S7 to S16 in section S5 of the online supplement plot the kernel density estimates of the parameters of the multivariate GARCH diffusion-driven volatility in mean model described above estimated by the MPM sampler (S = 100 and N = 250) using a 25% trimmed mean; the vertical lines show the true parameter values. The figures show that the model parameters are accurately estimated. Figure S6 in section S5 of the online supplement shows the inefficiency factors ($\widehat{\text{IF}}$) of the parameters over 10 simulated datasets with d = 20 dimensions and T = 100 time periods. The results suggest that the MPM sampler with only N = 250 particles and S = 100 particle filters is efficient because the $\widehat{\text{IF}}$ values of the parameters are reasonably small.

We now apply the methods to a sample of daily U.S. industry stock returns data obtained from the Kenneth French website,¹ using a sample from January 4th, 2021 to the 26th of May, 2021, a total of 100 observations. Table S7 in section S5 of the online supplement reports the $\widehat{\text{IF}}$, $\widehat{\text{TNIF}}$, and $\widehat{\text{RTNIF}}$ for the parameters of the 20 dimensional multivariate GARCH diffusion-driven volatility in mean model estimated using (a) the correlated PMMH, (b) the MPM sampler with the 25% trimmed mean of the likelihood estimate, and (c) the MPM sampler with the 50% trimmed mean of the likelihood estimate. The table shows that: (a) the MCMC chain obtained using the correlated PMMH with $N=100,000$ still gets stuck; (b) the CPU time for the correlated PMMH with $N=100,000$ particles is 26 times slower than the MPM algorithm with S = 250 and N = 100; (c) the MPM sampler with N = 250 and S = 100 is less efficient than the MPM sampler with N = 100 and S = 250. The MPM sampler with N = 100 and S = 250 can be made more efficient by running it on a high-performance computer cluster with many cores; (4) the best sampler in terms of ${\widehat{\text{TNIF}}}_{\text{MEAN}}$ and ${\widehat{\text{TNIF}}}_{\mathrm{\text{MAX}}}$ for estimating the 20 dimensional multivariate GARCH diffusion-driven volatility in mean model is the MPM sampler with a 50% trimmed mean with S = 250 and N = 100.

Table 4 reports the $\widehat{\text{IF}},\widehat{\text{TNIF}}$, and $\widehat{\text{RTNIF}}$ for the parameters of the d = 30 dimensional multivariate GARCH diffusion-driven volatility in mean model estimated using: (a) the MPM sampler with the 25% trimmed mean of the likelihood, and (b) the MPM sampler with 50% trimmed mean of the likelihood. We do not consider the correlated PMMH as it is very expensive to run for this high dimensional model. The table shows that: (a) increasing the number of particle filters from S = 100 to S = 200 reduces the $\widehat{\text{IF}}$, but increases the computational time substantially. However, the computational time can be reduced by using a high performance computer cluster with a larger number of cores. (b) the MPM sampler with N = 500 and S = 100 is less efficient than the MPM sampler with N = 250 and S = 200. If there is access to large computing resources, then we recommend increasing the number of particle filters S while keeping the number of particles N in each particle filter reasonably small.

3.4 DSGE Models

This section evaluates the effectiveness of the MPM samplers for estimating two nonlinear DSGE models: a stylized small-scale and a medium-scale DSGE model described in sections S8 and S9 of the online supplement. We use the disturbance filter method introduced in section 2.3 because the dimension of the state vector is much larger than that of the disturbances in these models. The state transition densities of the two nonlinear DSGE models are intractable. Section S3 of the online supplement describes the state space representation of these models, which are solved using Dynare.

3.4.1 Small Scale Nonlinear DSGE Model

This section reports on how well a number of PMMH samplers of interest estimate the parameters of a nonlinear (second order) small-scale DSGE model. This small-scale DSGE model is stylized and similar to the canonical three-equation New Keynesian model. The modeling framework follows Herbst and Schorfheide (2019) with nominal rigidities through price adjustment costs following Rotemberg (1982), and three structural shocks. There are four state variables, four control variables (8 variables combining the state and control variables) and 3 exogenous shocks in the model. We estimate this model using 3 observables, which means that the dimensions of y_t , the state vector, and the disturbance vector are 3, 8, and 3, respectively. Section S8 of the online supplement describes the model.

In this section, we use the delayed acceptance version of the MPM algorithm, calling it “delayed acceptance multiple PMMH” (DA-MPM), to speed up the computation. The delayed acceptance sampler (Christen and Fox 2005) avoids computing the expensive likelihood estimate if it is unlikely that the proposed draw will ultimately be accepted. The first accept-reject stage uses a cheap (or deterministic) approximation to the likelihood instead of the expensive likelihood estimate in the MH acceptance ratio. The particle filter is then used to estimate the likelihood only for a proposal that is accepted in the first stage; a second accept-reject stage ensures that detailed balance is satisfied with respect to the true posterior. We use the likelihood obtained from the central difference Kalman filter (CDKF) proposed by Norgaard, Poulsen, and Ravn (2000) in the first accept-reject stage of the delayed acceptance scheme. See section S2 of the online supplement for further details. The dimension of states is usually much larger than the dimension of the disturbances for the DSGE models. Two different sorting methods are also considered. The first sorts the state particles and the second sorts the disturbance particles using the sorting algorithm in Section 2.4.

The PMMH samplers considered are: (a) the MPM (0% trimmed mean) with the auxiliary disturbance particle filter (ADPF), disturbance sorting, and the correlation between random numbers used in estimating the log of estimated likelihoods set to ${\rho }_u=0.9$; (b) the MPM (10% trimmed mean) with the ADPF, disturbance sorting, and the ${\rho }_u=0.9$; (c) the MPM (25% trimmed mean) with the ADPF, disturbance sorting, and the ${\rho }_u=0.9$; (d) the MPM (25% trimmed mean) with the ADPF, disturbance sorting, and ${\rho }_u=0$, (e) the MPM (25% trimmed mean) with the ADPF, state sorting, and ${\rho }_u=0.9$; (f) the delayed acceptance MPM (DA-MPM) (25% trimmed mean) with ADPF, disturbance sorting, and the ${\rho }_u=0.9$; (g) the MPM (50% trimmed mean) with the ADPF, disturbance sorting, and the ${\rho }_u=0.9$; (h) the MPM (0% trimmed mean) with the bootstrap particle filter, disturbance sorting, and the ${\rho }_u=0.9$; (i) the correlated PMMH with ${\rho }_u=0.9$, the bootstrap particle filter, and disturbance sorting. Each sampler ran for $25,000$ iterations, with the initial 5000 iterations discarded as burn-in. We use the adaptive random walk proposal of Roberts and Rosenthal (2009) for $q({\theta }^{\prime }|\theta )$ and the adaptive scaling approach of Garthwaite, Fan, and Sisson (2015) to tune the Metropolis-Hastings acceptance probability to 20%. Section S8 of the online supplement gives details on model specifications and model parameters. We include measurement error in the observation equation. The measurement error variances are estimated together with other parameters. All the samplers ran on a high performance computer cluster with 20 CPU cores. The real dataset is obtained from Herbst and Schorfheide (2019), using data from 1983Q1 to 2013Q4, which includes the Great Recession with a total of 124 observations for each series.

Table 5 reports the relative time normalized inefficiency factor ($\widehat{\text{RTNIF}}$) of a PMMH sampler relative to the DA-MPM (25% trimmed mean) with the ADPF, disturbance sorting, and ${\rho }_u=0.9$ for the nonlinear small scale models with T = 124 observations. The computing time reported in the table is the time to run a single particle filter for the CPM and S particle filters for the MPM approach. The table shows that: (a) The MPM sampler (0% trimmed mean) with the ADPF, disturbance sorting, ${\rho }_u=0.9$, and N = 250 particles is more efficient than the MPM sampler (0% trimmed mean) with the bootstrap particle filter, disturbance sorting, ${\rho }_u=0.9$, and N = 250 particles with the MCMC chain obtained by the latter getting stuck. (b) The running time taken by the MPM method with S = 100 particles filters, each with N = 100 particles and disturbance sorting is 9.73 times faster than the CPM with $N=20,000$ particles. (c) In terms of ${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$, the performance of MPM samplers with 0%, 10%, 25%, and 50% trimmed means, ADPF, disturbance sorting, and ${\rho }_u=0.9$ are 3.15, 90.77, 103.04, and 91.50 times more efficient respectively, than the correlated PMMH with $N=20,000$ particles. (d) In terms of ${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$, the MPM (25% trimmed mean) with disturbance sorting, and ${\rho }_u=0.9$ is 1.38 times more efficient than the MPM (25% trimmed mean) method with state sorting method and performs similarly to the MPM (25% trimmed mean) approach with ADPF, disturbance sorting, and ${\rho }_u=0$. (e) The delayed acceptance MPM (25% trimmed mean) is slightly more efficient than the MPM (25% trimmed mean) approach. The delayed acceptance algorithm is 5 times faster on average because the target Metropolis-Hastings acceptance probability is set to 20% using the Garthwaite, Fan, and Sisson (2015) approach, but it has higher maximum $\widehat{\text{IF}}$ value. Table S9 in section S6 of the online supplement gives the details. (f) It is possible to use the tempered particle filter (TPF) of Herbst and Schorfheide (2019) within the MPM algorithm. However, the TPF is computationally expensive because of the tempering iterations and the random walk Metropolis-Hastings mutation steps. In summary, the best sampler with the smallest ${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$ and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ is the delayed acceptance MPM (DA-MPM) (25% trimmed mean) with ADPF, disturbance sorting, and ${\rho }_u=0.9$.

3.4.2 Nonlinear Medium Scale DSGE Model

Section 3.4.1 evaluates the performance of different PMMH samplers in estimating the stylized three equation New Keynesian model. A vast majority of studies however require a medium-scale model with a higher dimensional state space. This section, therefore, applies the MPM samplers with 50% trimmed mean for estimating the parameters of the second-order medium-scale DSGE model of Gust et al. (2017) discussed in section S9 of the online supplement.

The medium-scale DSGE model considered is similar to the models in Christiano, Eichenbaum, and Evans (2005) and Smets and Wouters (2007) with nominal rigidities in both wages and prices, habit formation in consumption, investment adjustment costs, and five structural shocks. We solve the model using perturbation methods and consider a second-order accurate solution. In the DSGE model examples, we consider, as before, the source of nonlinearity stems from the second-order interaction terms between endogenous and exogenous variables in the model. Our article does not consider occasionally binding constraints such as the zero lower bound (ZLB) on the nominal interest and solve the model using perturbation methods while focusing on the efficiency of the estimation algorithm. We restrict our sample to the period before the zero lower bound constraint on the nominal interest rate binds.

The baseline model in Gust et al. (2017) consists of 12 state variables, 18 control variables (30 variables combining the state and control variables) and 5 exogenous shocks. We estimate this model using five observables; that means that the dimensions of y_t , the state vector and the disturbance vector are 5, 30, and 5, respectively. We include measurement error in the observation equation and consider two cases. The measurement error variances are: (a) estimated and (b) fixed at the 25% of the sample variance of the observables over the estimation period. We use quarterly U.S. data from 1983Q1 to 2007Q4, with a total of 100 observations for each of the five series.

The PMMH samplers considered are: (a) the MPM (50% trimmed mean) with the auxiliary disturbance particle filter (ADPF), disturbance sorting, and the correlation between random numbers used in estimating the log of estimated likelihoods set to ${\rho }_u=0.9$; (b) the correlated PMMH of Deligiannidis, Doucet, and Pitt (2018). Each sampler ran for 60,000 iterations, with the initial 30,000 iterations discarded as burn-in. We use the adaptive random walk proposal of Roberts and Rosenthal (2009) for $q({\theta }^{\prime }|\theta )$. All the samplers ran on a high performance computer cluster with 20 CPU cores. Table 6 shows the $\widehat{\text{RTNIF}}$ of the correlated PMMH relative to the MPM sampler with 50% trimmed mean (ADPF, disturbance sorting, and ${\rho }_u=0.9$) for the nonlinear medium scale DSGE model for the cases where the measurement error variances are estimated and fixed at 25% sample variance of the observables. The computing time reported in the table is the time to run a single particle filter for the CPM and S = 250 particle filters for the MPM approach. The table shows that in terms of ${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$, the MPM sampler is 26 and 14 times more efficient than the correlated PMMH for the cases where the measurement error variances are estimated and fixed at 25% sample variance of the observables, respectively. Table S12 in section S7 of the online supplement gives the details on the inefficiency factors for each parameter.

Table 6 Comparing the performance of different PMMH samplers using different numbers of particle filters S and different numbers of particles N in each particle filter for estimating the medium scale DSGE model using the U.S. quarterly dataset from 1983Q1 to 2007Q4.

	I	II	I actual	I	II	I actual
${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}}$	1.00	2.74	796.82	1.00	1.48	383.31
${\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}}$	1.00	26.51	1338.66	1.00	14.37	643.96
${\widehat{\text{IF}}}_{\text{MEAN}}$	1.00	2.78	428.71	1.00	0.96	275.34
${\widehat{\text{RTNIF}}}_{\text{MEAN}}$	1.00	26.92	720.22	1.00	9.31	462.57
Time	1.00	9.68	1.68	1.00	9.68	1.68

NOTE: Sampler I: MPM (50% trimmed mean, ${\rho }_u=0.9$, disturbance sorting, ADPF) with S = 250 and N = 100. Sampler II: Correlated PMMH with $N=50,000$. The first three columns are for the case when the measurement error variances are estimated. The ${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{IF}}}_{\text{MEAN}}$, and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ entries in columns headed I–II are relative to the entries in column I. The entries in column “I actual” are the actual ${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{IF}}}_{\text{MEAN}}$, and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ for estimator I. The last three columns are for the case when the measurement error variances are fixed at 25% of the sample variance of the observables. The ${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{IF}}}_{\text{MEAN}}$, and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ entries in columns headed I to II are relative to the entries in column I. The entries in column “I actual” are the actual ${\widehat{\text{IF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{RTNIF}}}_{\mathrm{\text{MAX}}},{\widehat{\text{IF}}}_{\text{MEAN}}$, and ${\widehat{\text{RTNIF}}}_{\text{MEAN}}$ for estimator I.

Table S11 in section S7 of the online supplement gives the posterior means, 2.5%, and 97.5% quantile estimates of all the parameters in the medium scale DSGE model estimated using the MPM sampler with 50% trimmed mean (${\rho }_u=0.9$, disturbance sorting, ADPF). The measurement error variances are fixed to 25% of the sample variance of the observables. The table shows that the standard errors are relatively small for all parameters indicating that the parameters are estimated accurately. Figure 1 shows the kernel density estimates of some of the parameters of the medium scale DSGE model estimated using the MPM with 50% trimmed mean (ADPF, disturbance sorting, and ${\rho }_u=0.9$) and the correlated PMMH. The figure shows that the approximate posteriors obtained using the MPM with 50% trimmed mean is very close to the exact posteriors obtained using the correlated PMMH.

Fig. 1 Kernel density estimates of the posterior densities of some parameters of the medium scale DSGE model for the U.S. quarterly dataset from 1983Q1 to 2007Q4 estimated using: (a) MPM with 50% trimmed mean (${\rho }_u=0.9$, disturbance sorting, ADPF) with S = 250 and N = 100; (b) correlated PMMH with $N=50,000$. The measurement error variances are fixed to 25% of the variance of the observables.

We also consider a variation of the medium scale DSGE model of Gust et al. (2017) described in section S9 of the online supplement. To study how well our approach performs with increasing state dimension, we estimate the model with the canonical definition of the output gap defined as the difference between output in the model with nominal rigidities and output in the model with flexible prices and wages. This variation is identical to the baseline model in Gust et al. (2017) except for the Taylor rule specification, which now uses the canonical definition of the output gap. This extension of the medium scale DSGE model consists of 21 state variables, 30 control variables (51 variables combining the state and control variables) and 5 exogenous shocks. We estimate this model using five observables: this means that the dimensions of y_t , the state vector and the disturbance vector are 5, 51, and 5, respectively. The measurement error variances are fixed at 25% of the variance of the observables over the estimation period. Table S13 of the online supplement shows the posterior means, 2.5% and 97.5% quantiles of the parameters of the extended medium scale DSGE model estimated using the MPM algorithm with a 50% trimmed mean (ADPF, disturbance sorting, and ${\rho }_u=0.9$) with S = 250 and N = 100. The table shows that the $\widehat{\text{IF}}$s of the parameters are similar to the simpler medium scale DSGE model, indicating the performance of the MPM algorithm does not deteriorate with the extended model having higher state dimensions.

4 Summary and Conclusions

The article proposes a general pseudo marginal approach (MPM) for estimating the posterior density of the parameters of complex and high-dimensional state-space models. It is especially useful when the single bootstrap filter is inefficient, while the more efficient auxiliary particle filter cannot be used because the state transition density is computationally intractable. The MPM method is a general extension of the PMMH method and makes four innovations: (a) it is based on the mean or trimmed mean of unbiased likelihood estimators; (b) it uses a novel block version of PMMH that works with multiple particle filters; (c) an auxiliary disturbance particle filter sampler is proposed to estimate the likelihood which is especially useful when the dimension of the states is much larger than the dimension of the disturbances; (d) a fast Euclidean sorting algorithm is proposed to preserve the correlation between the logs of the estimated likelihoods at the current and proposed parameter values. The MPM samplers are applied to complex multivariate GARCH diffusion-driven volatility in mean and DSGE models and are shown to work well.

The source of nonlinearity we consider in estimating DSGE models is due to the interaction terms in the model solution. The article aims to demonstrate the proposed algorithm’s efficacy in estimating complicated state space models with a high-dimensional state space. There are other sources of nonlinearities in the DSGE model that we do not consider in the present article. For instance, the model does not consider the zero lower bound on nominal interest rates. With occasionally binding constraints, such frameworks require global methods to solve the model before evaluating the likelihood. However, solution methods using global techniques are time-intensive. To reduce the computational time, we focus exclusively on models which can be solved using perturbation methods. Given a solution technique, the proposed algorithm can be used to estimate nonlinear DSGE models incorporating important features such as occasionally binding constraints. The other feature we do not address is the case of big shocks such as the Covid-19 pandemic. Modeling rare-disaster shocks is challenging when solving models using perturbation methods. Fernández-Villaverde and Levintal (2018) show that even at a third order, the model solution may suffer from accuracy issues when considering big shocks. These are interesting avenues for future work.

Supplementary Materials

The article has an online supplement. The commented Matlab code and readme.txt files are available at https://github.com/davidgun1000/Block-Correlated-Pseudo-Marginal to enable the reader to replicate all the empirical examples in the article.

Acknowledgments

We would like to thank the anonymous associate editor and two referees for their comments which helped improve both the presentation and technical content of the article.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The research of Robert Kohn and David Gunawan was partially supported by an ARC Center of Excellence grant CE140100049.

Notes

1 http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/datalibrary.html