The Block-Correlated Pseudo Marginal Sampler for State Space Models

We use the colon notation for collections of variables, i.e., $a_{t}^{r:s}:=\left(a_{t}^{r},...,a_{t}^{s}\right)$ and for $t\leq u$, $a_{t:u}^{r:s}:=\left(a_{t}^{r:s},...,a_{u}^{r:s}\right)$. Suppose $\left\{\left(Z_{t},Y_{t}\right),t\geq 0\right\}$ is a stochastic process, with parameter $\theta$, 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\left(z_{0}|\theta\right)$ the density of $Z_{0},$ $p\left(z_{t}|z_{t-1},\theta\right)$ the density of $Z_{t}$ given $Z_{0:t-1}=z_{0:t-1}$ for $t\geq 1$, and $p\left(y_{t}|z_{t},\theta\right)$ is the density of $Y_{t}$ given $Z_{0:t}=z_{0:t}$, $Y_{1:t-1}=y_{1:t-1}$.

The objective of Bayesian inference is to obtain the posterior distribution of the model parameters $\theta$ and the latent states $z_{0:T}$, given the observations $y_{1:T}$ and a prior distribution $p\left(\theta\right)$; i.e.,

$\displaystyle p\left(\theta,z_{0:T}|y_{1:T}\right)$

$\displaystyle=$

$\displaystyle 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\left(y_{1:T}\right)=\int\int p\left(y_{1:T}|\theta,z_{0:T}\right)p\left(z_{0:T}|\theta\right)p\left(\theta\right)dz_{0:T}d\theta$$

is the marginal likelihood. The likelihood

$$p\left(y_{1:T}|\theta\right)=\int p\left(y_{1:T}|\theta,z_{0:T}\right)p\left(z_{0:T}|\theta\right)dz_{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 non-linear and non-Gaussian state space models and the non-linear DSGE models described in section 3.4.

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

The non-linear (second-order) DSGE models considered in section 3.4 lie within the class of general non-linear 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 Eq. (1) becomes

$$p\left(\theta,\epsilon_{1:T},z_{0}|y_{1:T}\right)\propto\prod_{t=1}^{T}p\left(y_{t}|\epsilon_{1:t},\theta,z_{0}\right)p\left(\epsilon_{t}\right)p\left(z_{0}|\theta\right)p\left(\theta\right),$$

(2)

which Murray et al., (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\left(z_{t-1},\epsilon_{t};\theta\right)\rightarrow z_{t}$. This gives us a state trajectory $z_{0:T}$ from any sample $\left(\epsilon_{1:T},z_{0}\right)$. In the disturbance state-space model the target becomes the posterior distribution over the parameters $\theta$, the latent noise variables $\epsilon_{1:T}$ and the initial state $z_{0}$, rather than $\left(\theta,z_{0:T}\right)$.

This section outlines the standard PMMH scheme which carries out MCMC on an expanded space using an unbiased estimate of the likelihood. Our paper focuses on MCMC based on the parameters, but it is straightforward to also use it to obtain posterior inference for the unobserved states; see, e.g., Andrieu et al., (2010). Let $u$ consist of all the random variables required to compute the unbiased likelihood estimate $\widehat{p}_{N}\left(y|\theta,u\right)$, with $p(u)$ the density of $u$; let $p\left(\theta\right)$ be the prior of $\theta$. The joint posterior density $\theta$ and $u$ is

$$p\left(\theta,u|y_{1:T}\right)=\widehat{p}_{N}\left(y_{1:T}|\theta,u\right)p\left(\theta\right)p\left(u\right)/p\left(y_{1:T}\right),$$

so that

$$p\left(\theta|y_{1:T}\right)=\int p\left(\theta,u|y_{1:T}\right)du=p\left(y_{1:T}|\theta\right)p\left(\theta\right)/p\left(y_{1:T}\right)$$

is the posterior of $\theta$ and $\intop\widehat{p}_{N}\left(y_{1:T}|\theta,u\right)p\left(u\right)du=p\left(y_{1:T}|\theta\right)$ because the likelihood estimate is unbiased. We can therefore sample from the posterior density $p\left(\theta|y_{1:T}\right)$ by sampling $\theta$ and $u$ from $p\left(\theta,u|y_{1:T}\right)$. The subscript $N$ indicates the number of particles used to estimate the likelihood.

Let $q\left(\theta^{{}^{\prime}}|\theta\right)$ be the proposal density for $\theta^{\prime}$ with the current value $\theta$ and $q\left(u^{{}^{\prime}}|u\right)$ the proposal density for $u^{{}^{\prime}}$ given $u$. We always assume that $q\left(u^{{}^{\prime}}|u\right)$ satisfies the reversibility condition

$$q\left(u^{{}^{\prime}}|u\right)p\left(u\right)=q\left(u|u^{{}^{\prime}}\right)p\left(u^{{}^{\prime}}\right);$$

it is clearly satisfied by the standard PMMH where $q\left(u^{{}^{\prime}}|u\right)=p\left(u^{{}^{\prime}}\right)$. We generate a proposal $\theta^{{}^{\prime}}$ from $q\left(\theta^{{}^{\prime}}|\theta\right)$ and $u^{{}^{\prime}}$ from $q\left(u^{{}^{\prime}}|u\right)$, and accept both proposals with probability

	$\displaystyle\alpha\left(\theta,u;\theta^{{}^{\prime}},u^{{}^{\prime}}\right)$	$\displaystyle=$	$\displaystyle\min\left(1,\frac{\widehat{p}_{N}\left(y|\theta^{{}^{\prime}},u^{{}^{\prime}}\right)p\left(\theta^{{}^{\prime}}\right)p\left(u^{{}^{\prime}}\right)q\left(\theta|\theta^{{}^{\prime}}\right)q\left(u|u^{{}^{\prime}}\right)}{\widehat{p}_{N}\left(y|\theta,u\right)p\left(\theta\right)p\left(u\right)q\left(\theta^{{}^{\prime}}|\theta\right)q\left(u^{{}^{\prime}}|u\right)}\right)$		(3)
		$\displaystyle=$	$\displaystyle\min\left(1,\frac{\widehat{p}_{N}\left(y|\theta^{{}^{\prime}},u^{{}^{\prime}}\right)p\left(\theta^{{}^{\prime}}\right)q\left(\theta|\theta^{{}^{\prime}}\right)}{\widehat{p}_{N}\left(y|\theta,u\right)p\left(\theta\right)q\left(\theta^{{}^{\prime}}|\theta\right)}\right).$

The expression in Eq. (3) is identical to a standard Metropolis-Hastings algorithm, except that estimates of the likelihood at the current and proposed parameters are used. Andrieu and Roberts, (2009) show that the resulting PMMH algorithm has the correct invariant distribution regardless of the variance of the estimated likelihoods. However, the performance of the PMMH approach crucially depends on the number of particles $N$ used to estimate the likelihood. The variance of the log of the estimated likelihood should be between 1 and 3 depending on the quality of the proposal for $\theta$; see Pitt et al., (2012) and Sherlock et al., (2015). In many applications of the non-linear state space models considered in this paper, it is computationally very expensive to ensure that the variance of the log of the estimated likelihood is within the required range. Section 2.3 discusses the new PMMH sampler, which we call the multiple PMMH algorithm (MPM), that builds on and extends, the CPM of Deligiannidis et al., (2018) and BPM of Tran et al., (2016).

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 block-correlated PMMH that works with multiple particle filters to induce a high correlation between successive logs of the estimated likelihoods is also discussed. We now define the joint target density of $\theta$ and $\widetilde{u}=\left(u_{1},...,u_{S}\right)$ as

$$p\left(\theta,\widetilde{u}|y_{1:T}\right)\propto\overline{\widehat{p}}_{N}\left(y|\theta,\widetilde{u}\right)p\left(\theta\right)\prod_{s=1}^{S}p\left(u_{s}\right),$$

(4)

where $\overline{\widehat{p}}_{N}\left(y|\theta,u\right)$ is the likelihood estimates obtained from $S$ particle filters discussed below. We update the parameters $\theta$ jointly with a randomly selected block $u_{s}$ in each MCMC iteration, with $\Pr\left(S=s\right)=1/S$ for any $s=1,...,S$. The selected block $u_{s}$ is updated using $u_{s}^{{}^{\prime}}=\rho_{u}u_{s}+\sqrt{1-\rho_{u}^{2}}\eta_{u}$, where $\rho_{u}$ is the non-negative correlation between the random numbers $u_{s}$ and $u_{s}^{{}^{\prime}}$ and $\eta_{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, e.g. Johnson et al., (2013). Using this scheme, the acceptance probability is

$$\alpha\left(\theta,\widetilde{u};\theta^{{}^{\prime}},\widetilde{u}^{{}^{\prime}}\right)=\min\left(1,\frac{\overline{\widehat{p}}_{N}\left(y|\theta^{{}^{\prime}},\widetilde{u}^{{}^{\prime}}=\left(u_{1},...,u_{s-1},u_{s}^{{}^{\prime}},u_{s+1},...,u_{S}\right)\right)p\left(\theta^{{}^{\prime}}\right)q\left(\theta|\theta^{{}^{\prime}}\right)}{\overline{\widehat{p}}_{N}\left(y|\theta,\widetilde{u}=\left(u_{1},...,u_{s-1},u_{s},u_{s+1},...,u_{S}\right)\right)p\left(\theta\right)q\left(\theta^{{}^{\prime}}|\theta\right)}\right).$$

(5)

We now discuss approaches to obtain a likelihood estimate from multiple particle filters which we then use in the MPM algorithm. Suppose that $S$ particle filters are run in parallel. Let $\widehat{p}_{N}\left(y|\theta,u_{s}\right)$ be the unbiased estimate of the likelihood obtained from the $s$th particle filter, for $s=1,...,S$. The first approach takes the average of the $S$ unbiased likelihood estimates from the particle filter. The resulting likelihood estimate

$$\overline{\widehat{p}}_{N}\left(y|\theta,u\right)=\sum_{s=1}^{S}\widehat{p}_{N}\left(y|\theta,u_{s}\right)/S$$

(6)

is also unbiased (Sherlock et al., 2017b, ). The second approach is to take the trimmed mean as the likelihood estimate in MPM. 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 approximates the exact posterior well.

It is important that the logs of the likelihood estimates evaluated at the current and proposed values of $\theta$ and $\widetilde{u}$ are highly correlated to reduce the variance of $\log\overline{\widehat{p}}_{N}\left(y|\theta^{{}^{\prime}},\widetilde{u}^{{}^{\prime}}\right)-\log\overline{\widehat{p}}_{N}\left(y|\theta,\widetilde{u}\right)$ 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 $\theta$ 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 et al., , 2018). However, such simple sorting is unavailable for multidimensional state particles.

This section discusses Algorithm 2, 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 S1 in section S1 of the online supplement. The algorithm is written 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 normalised weights $\overline{w}_{t}^{1:N}$ as the inputs; it outputs the sorted particles $\widetilde{z}_{t}^{1:N}$, sorted weights $\widetilde{\overline{w}}_{t}^{1:N}$, and sorted indices $\zeta_{1:N}$. Let $z_{t}^{i}$ be $d$-dimensional particles at a time $t$ for the $i$th particle, $z_{t}^{i}=\left(z_{t,1}^{i},...,z_{t,d}^{i}\right)^{\top}$. Let $d\left(z_{t}^{j},z_{t}^{i}\right)$ be the Euclidean distance between two multidimensional particles. The first step is to calculate the means of the multidimensional particles $\overline{z}_{t}^{i}=(\nicefrac{{1}}{{d}})\sum_{k=1}^{d}z_{t,k}^{i}$ at time $t$ for $i=1,...,N$. The first sorting index, $\zeta_{1}$, is the index of the multidimensional particle having the smallest value of $\overline{z}_{t}^{i}$ for $i=1,...,N$. Therefore, the first sorted particle $\widetilde{z}_{t}^{1}$ is $z_{t}^{\zeta_{1}}$ with its associated weight $\widetilde{\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\left(\widetilde{z}_{t}^{1},z_{t}^{i}\right)$, for all $i\neq\zeta_{1}$. The next step is to sort the particles and weights according to the Euclidean distance from smallest to largest to obtain the sorted particles $\widetilde{z}_{t}^{2:N}$, sorted weights $\widetilde{\overline{w}}_{t}^{2:N}$, and sorted indices $\zeta_{2:N}$.

Algorithm 2 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, $\zeta_{2}$, we need to calculate the Euclidean distance between the first selected particle and the remaining multidimensional particles. Similarly, to select the third sorted index, $\zeta_{3}$, 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 long time series. In contrast, Algorithm 2 only needs to calculate the Euclidean distance of the first selected particle and the remaining particles once. Deligiannidis et al., (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 2. Table 1 shows the CPU time (in seconds) of the three sorting methods for different numbers of particles $N$ and state dimensions $d$. The table shows that 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 2 is $290$ and $335$ times faster than the Choppala et al., (2016) and the Hilbert sorting methods, respectively

This section discusses the auxiliary disturbance particle filter we use to obtain the estimates of the likelihood in the MPM sampler described in section 2.3. 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}=\Phi\left(z_{0},\epsilon_{1};\theta\right)$, where $z_{0}$ is the initial state vector with density $p(z_{0}|\theta)$, and $z_{t}=F\left(z_{t-1},\epsilon_{t};\theta\right)(t\geq 2)$, where $\epsilon_{t}$ is a $n_{e}\times 1$ vector of normally distributed latent noise with density $p\left(\epsilon_{t}\right).$ Murray et al., (2013) express the standard state-space model in terms of the latent noise variables $\epsilon_{1:T}$, and call

$$\epsilon_{t}\sim p\left(\epsilon_{t}\right),y_{t}|\epsilon_{1:t}\sim p\left(y_{t}|\epsilon_{1:t},z_{0};\theta\right)=p\left(y_{t}|z_{t};\theta\right),\,\,t=1,...,T,$$

the disturbance state-space model. We note that the conditional distribution of $y_{t}$ depends on all the latent error variables, $\epsilon_{1:t}.$

We now discuss the proposal for $\epsilon_{t}$ used in the disturbance filter. The defensive mixture proposal density (Hesterberg, , 1995) is

$$m\left(\epsilon_{t}|u_{\epsilon,t},\theta\right)=\pi p\left(\epsilon_{t}|\theta\right)+(1-\pi)q\left(\epsilon_{t}|\theta,y_{1:T}\right),\quad{\rm with}\quad 0<\pi\ll 1,$$

(7)

where $u_{\epsilon,t}$ is the vector random variable used to generate the particles $\epsilon_{t}$ given $\theta$. If the observation density $p\left(y_{t}|z_{t},\theta\right)$ is bounded, then Eq. (7) guarantees that the weights are bounded in the disturbance particle filter algorithm defined in Eq. (S2) of section S1 of the online supplement. In practice, we take $q\left(\epsilon_{t}|\theta,y_{1:T}\right)=N\left(\epsilon_{t}|\widehat{\mu}_{t},\widehat{\Sigma}_{t}\right)$ and set $\pi=0.05$ in all the examples in section 3; see algorithm 4. If $\pi=1$ in Eq. (7), then $m\left(\epsilon_{t}\right)=p\left(\epsilon_{t}\right)$ 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 4 discusses how we obtain $\widehat{\mu}_{t}$ and the covariance matrix $\widehat{\Sigma}_{t}$, for $t=1,...,T$, at each iteration of the MPM algorithm. It is computationally cheap to obtain $\widehat{\mu}_{t}$ and $\widehat{\Sigma}_{t}$ for $t=1,...,T$ because they are obtained 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,...,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 initialise the mean $\widehat{\mu}_{t}$ and the covariance matrix $\widehat{\Sigma}_{t}$, for $t=1,...,T$. Algorithm 3 gives the MPM algorithm with ADPF.

We use the inefficiency factor (IF) (also called the integrated autocorrelation time)

$$\textrm{IF}_{\psi}:=1+2\sum_{j=1}^{\infty}\rho_{\psi}\left(j\right),$$

to measure the inefficiency of a PMMH sampler at estimating the posterior expectation of a univariate function $\psi\left(\theta\right)$ of $\theta$; here, $\rho_{\psi}\left(j\right)$ is the $j$th autocorrelation of the iterates $\psi\left(\theta\right)$ in the MCMC chain after it has converged to its stationary distribution. We estimate the $\textrm{IF}_{\psi}$ using the CODA R package of Plummer et al., (2006). A low value of the $\textrm{IF}_{\psi}$ 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 $\theta$ based on $\textrm{IF}_{\psi}$ is the time normalised inefficiency factor (TNIF) defined as $\textrm{TNIF}_{\psi}:=\textrm{IF}_{\psi}\times\textrm{CT}.$ For a given sampler, let $\textrm{IF}_{\psi,\textrm{MAX}}$ and $\textrm{IF}_{\psi,\textrm{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.

This section investigates (1) the efficiency of different approaches for obtaining likelihood estimates from $S$ particle filters, and (2) the performance of different PMMH samplers for estimating a single parameter $\theta$, regardless of the dimension of the states. We consider the linear Gaussian state space model discussed in Guarniero et al., (2017) and Deligiannidis et al., (2018), where $\left\{X_{t};t\geq 1\right\}$ and $\left\{Y_{t};t\geq 1\right\}$ are ${\cal R}^{d}$ valued with

	$\displaystyle Y_{t}$	$\displaystyle=$	$\displaystyle X_{t}+W_{t},$
	$\displaystyle X_{t+1}$	$\displaystyle=$	$\displaystyle A_{\theta}X_{t}+V_{t+1},$

with $X_{1}\sim{\cal N}\left(0_{d},I_{d}\right)$, $V_{t}\sim N\left(0_{d},I_{d}\right)$, $W_{t}\sim N\left(0_{d},I_{d}\right)$, and $A_{\theta}^{i,j}=\theta^{|i-j|+1}$ for $i<j$; the true value of $\theta$ 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 shows the 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$. (2) 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.

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 $\widetilde{u}$, where $\widetilde{u}=\left(u_{1},...,u_{S}\right)$ and $\widetilde{u}^{{}^{\prime}}$ and obtained $\log\overline{\widehat{p}}_{N}\left(y|\theta,\widetilde{u}^{\left(s\right)}\right)$ and $\log\overline{\widehat{p}}_{N}\left(y|\theta^{{}^{\prime}},\widetilde{u}^{\left(s\right)^{\prime}}\right)$ for the MPM approaches and computed their sample correlations.

Figure S2 in Section S4 of the online supplement reports the correlation estimates of the log of estimated likelihood obtained using different MPM approaches. The figure show that: (1) 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; (2) when the current parameter $\theta$ 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 different sampling schemes for estimating the parameter $\theta$ of the linear Gaussian state space model. In all examples, we run the samplers for $25000$ iterations, with the initial $5000$ iterations discarded as warm up. The particle filter and the parameter samplers are implemented in Matlab running on $20$ CPU cores, a high performance computer cluster. 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 2 shows the trace plots of the parameter $\theta$ 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: (1) 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$; (2) the 10% trimmed mean approach requires at least $S=100$ and $N=250$ to make the MCMC chain mix well; (3) 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 1 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{\textrm{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). This suggests the usefulness of the MPM (blocking) sampling scheme.

Table 3 shows the $\widehat{\textrm{IF}}$, $\widehat{\textrm{TNIF}}$, and $\widehat{\textrm{RTNIF}}$ values for the parameter $\theta$ in the linear Gaussian state space model with $d=10$ dimensions and $T=300$ time periods estimated using the following five different MCMC samplers: (1) the correlated PMMH of Deligiannidis et al., (2018), (2) the MPM with 5% trimmed mean, (3) the MPM with 10% trimmed mean, (4) the MPM with 25% trimmed mean, and (5) the MPM with 50% trimmed mean. In this paper, 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 points. (1) The correlated PMMH requires more than $50000$ particles to improve the mixing of the MCMC chain for the parameter $\theta$. (2) The CPU time for running a single particle filter with $N=50000$ 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. (3) The MPM allows us to use much smaller number of particles for each independent PF and these multiple PFs can be run independently. (4) The $\widehat{\textrm{IF}}$ values for the parameter $\theta$ estimated using the correlated PMMH with $N=50000$ particles is $11.05$, $41.93$, $46.12$, and $56.37$ times larger than the 5%, 10%, 25%, and 50% trimmed means approaches with $S=100$ and $N=500$ particles, respectively. (5) In terms of $\widehat{\textrm{RTNIF}}$, the 5%, 10%, 25%, and 50% trimmed means with $S=100$ and $N=500$ are $45.27$, $171.52$, $188.67$, and $230.60$ times smaller than the correlated PMMH with $N=50000$ particles. (6) The best sampler for this example is the 50% trimmed mean approach with $S=100$ and $N=250$ particles. Table LABEL:tabledim10T200 in section S4 of the online supplement gives similar results for the 10 dimensional linear Gaussian state space model with $T=200$. Figure 3 shows the kernel density estimates of the parameter $\theta$ 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 likelihood. The figure shows 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 S1 in Section S4 of the online supplement. Figure S1 also shows that the approximate posterior obtained using the MPM with trimmed mean of the likelihood is very close to the exact posterior obtained using the correlated PMMH.

In summary, the example suggests that: (1) 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 of the likelihood reduce the variance substantially when $S$ and/or $N$ increases; (2) the 25% and 50% trimmed means approaches give the smallest variance of the log of the estimated likelihood for all cases; (3) the MPM method with 50% trimmed mean is best at maintaining the correlation between the logs of the estimated likelihoods in successive iterates; (4) the MPM (blocking) method is more efficient than the MPM (no blocking) for the same values of $S$ and $N$; (5) the approximate approaches that use the trimmed means of the likelihood gives accurate approximations to the true posterior; (6) 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; (7) 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.

This section investigates the performance of the proposed MPM samplers for estimating large multivariate stochastic volatility in mean models (see, e.g., Carriero et al., , 2018; Cross et al., , 2021), where each of the log volatility processes follows a GARCH (Generalized Autoregressive Conditional Heteroskedasticity) diffusion process (Shephard et al., , 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}=\log P_{t}-\log P_{t-1}$ as the log-return of the stocks. Let $h_{i,t}$ be the log-volatility process of the $i$th stock at time $t$. We also define, $h_{\cdotp,t}=\left(h_{1,t},...,h_{d,t}\right)^{\top}$ and $h_{i,\cdotp}=\left(h_{i,1},...,h_{i,T}\right)^{\top}$. The model for $y_{t}$ is

$$y_{t}=Ah_{\cdot,t}+\epsilon_{t},\;\epsilon_{t}\sim N\left(0,\Sigma_{t}\right),$$

(8)

where $A$ is a $d\times d$ matrix that captures the effects of log-volatilities on the levels of the variables. The time-varying error covariance matrix $\Sigma_{t}$ depends on the unobserved latent variables $h_{\cdotp,t}$ such that

$$\Sigma_{t}:=\textrm{diag}\left(\exp\left(h_{1,t}\right),...,\exp\left(h_{d,t}\right)\right).$$

(9)

Each log-volatility $h_{i,t}$ is assumed to follow a continuous time GARCH diffusion process (Shephard et al., , 2004) satisfying

$$dh_{i,t}=\left\{\alpha\left(\mu-\exp\left(h_{i,t}\right)\right)\exp\left(-h_{i,t}\right)-\frac{\tau^{2}}{2}\right\}dt+\tau dW_{i,t},$$

(10)

for $i=1,...,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 (10) by placing $M-1$ evenly spaced points between times $t$ and $t+1$. We denote the intermediate volatility components by $h_{i,t,1},...,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, for example, Stramer and Bognar, (2011), pg. 234)

$$h_{i,t,j+1}|h_{i,t,j}\sim N\left(h_{i,t,j}+\left\{\alpha\left(\mu-\exp\left(h_{i,t,j}\right)\right)\exp\left(-h_{i,t,j}\right)-\frac{\tau^{2}}{2}\right\}\delta,\tau^{2}\delta\right),$$

(11)

for $j=0,...,M-1$, where $\delta=1/M$. The initial state of $h_{i,t}$ is assumed normally distributed $N(0,1)$ for $i=1,...,d$. For this example, we assume for simplicity that the parameters $\mu$, $\alpha$, 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.