{'fraction_non_alphanumeric': 0.035522420787194385, 'fraction_numerical': 0.016582611555750467, 'mean_word_length': 4.316188255864782, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Inference for State Space Models using Block and Correlated Pseudo Marginal Methods

This article addresses the problem of efficient Bayesian inference in dynamic systems using particle methods and makes a number of contributions. First, we develop a correlated pseudo-marginal (CPM) approach for Bayesian inference in state space (SS) models that is based on filtering the disturbances, rather than the states. This approach is useful when the state transition density is intractable or inefficient to compute, and also when the dimension of the disturbance is lower than the dimension of the state. Second, we propose a block pseudo-marginal (BPM) method that uses as the estimate of the likelihood the average of G independent unbiased estimates of the likelihood. We associate a set of underlying uniform of standard normal random numbers used to construct each of the individual unbiased likelihood estimates and then use component-wise Markov Chain Monte Carlo to update the parameter vector jointly with one set of these random numbers at a time. This induces a correlation of approximately 1-1/G between the logs of the estimated likelihood at the proposed and current values of the model parameters. Third, we show for some non-stationary state space models that the BPM approach is much more efficient than the CPM approach, because it is difficult to translate the high correlation in the underlying random numbers to high correlation between the logs of the likelihood estimates. Although our focus has been on applying the BPM method to state space models, our results and approach can be used in a wide range of applications of the PM method, such as panel data models, subsampling problems and approximate Bayesian computation.


Introduction
It is challenging to carry out Bayesian inference for the model parameters in SS models as the likelihood is often intractable. Because this likelihood at any given value of the parameters can be estimated unbiasedly using a particle filter (PF) (Del Moral, 2004), Andrieu et al. (2010) propose using the pseudo-marginal (PM) Markov chain Monte Carlo (MCMC) based on an unbiased estimate of the likelihood to sample from the posterior distribution of the unknown model parameters. Pitt et al. (2012) show that to obtain a good tradeoff between computational complexity and MCMC mixing, the number of particles used in the PF should be such that the variance of the log of the estimated likelihood (which we call σ 2 ) is around one and that the inefficiency of the PM scheme increases exponentially with σ 2 . Furthermore, the variance of the log of the estimated likelihood increases linearly with the sample size T and is inversely proportional to the number of particles for a given sample size. Hence, the number of particles required to keep the variance around 1 is O(T 2 ). Dahlin et al. (2015a) and Deligiannidis et al. (2016) recently propose the CPM algorithm which correlates the log of the estimated likelihood at the proposed and current values of the model parameters by correlating the underlying standard normal random numbers used to construct the estimates of the likelihood. Introducing this correlation into the likelihood estimates reduces the variance of the difference in the logs of the estimated likelihoods which appears in the Metropolis-Hastings acceptance ratio. Deligiannidis et al. (2016) show that the number of particles required by CPM in each MCMC iteration is O(T 3/2 ). Tran et al. (2016) propose an alternative PM approach, called the BPM approach, which divides the uniform or standard normal random numbers into blocks and then updates the unknown parameters jointly with one block of these random numbers, in each MCMC iteration. Tran et al. (2016) show that the number of particles required by BPM is O(T 3/2 ) if the likelihood is estimated using Monte Carlo, and is O(T 7/6 ) if the likelihood is estimated based on randomized quasi-Monte Carlo.
In some applications of SS models such as marine biogeochemical and economic models (Murray et al., 2013;Hall et al., 2014), it is difficult to estimate the likelihood unbiasedly using the standard particle filters that filter the state, because the state transition density is intractable or computationally expensive to compute. Murray et al. (2013) express this intractable transition density using the disturbances in the state transition equation, and reformulate the SS model by using the disturbances as the new state, who transition density is now tractable. It is possible, and sometimes more efficient, to carry out a particle filter that is based on filtering these disturbances. For example, in some SS models such as the Stochastic Volatility models in Section 6.1, the disturbance dimension is lower than the state dimension. We follow Murray et al. (2013) and refer to the SS model based on the disturbances as the Disturbance SS Model and the PF based on disturbances as the disturbance PF.
The first contribution of this article is to extend the CPM for use with the disturbance PF. The CPM algorithm operates on the space of states in such a way that the correlation in the likelihood values is preserved by maintaining similarity in the particle states. Our article proposes maintaining similarity in the disturbance space rather than in the state space. The motivating rationale is that if the disturbances are close to one another then the corresponding particle states will also be close to one another and hence the likelihood correlation will be preserved. Inducing strong positive correlation in the log of the estimated likelihood values leads to non-sticky mixing of the MCMC chain with fewer particles. In addition, it is more convenient and efficient to work with disturbances when their dimension is lower than that of the original states.
The CPM method involves a sorting procedure that can be computationally expensive; see Section 4. The second contribution of this article is to propose a highly efficient BPM scheme for Bayesian inference in the general class of time series SS models. The proposed BPM method extends the main idea in Tran et al. (2016) to time series models for which the unbiased likelihood estimator is an average of likelihood estimates obtained by multiple independent PFs. These PF's are run in parallel on multiple processors, with the unknown model parameters and one block of the uniform or standard normal random numbers used in one of the PFs updated jointly. Using the power of multiple-processor architecture in the PM context has been explored recently by Drovandi (2014), but he did not incorporate the blocking idea. This extension of the BPM approach in Tran et al. (2016) leads to a surprisingly efficient Bayesian inference approach for SS models. The theory in Tran et al.
(2016) makes it possible to obtain a desired correlation for the log of the estimated likelihoods to a sufficient accuracy, unlike the CPM method, and depends only on the number of independent PFs used.
Although our focus has been on applying the BPM method to SS models, our results and approach apply equally to a number of other applications of the PM method such as panel data models, subsampling problems and approximate Bayesian computation.
The rest of the paper is organised as follows: Section 2 sets out the notation for the conventional SS model. Section 3 presents the disturbance PF methodology.
Section 4 considers the proposed CPM method for disturbance particle filtering. 4 Section 5 presents the BPM method. Section 6 illustrates the methodology using simulated examples, and section 7 applies the methodology to real data. Section 8 concludes. There are two appendices. The first outlines how to estimate the gradient and hessian of the likelihood, and the second gives further details on the analysis of the simulated and real examples.

Conventional SS model
The SS model, shown in Fig. 1, describes the evolution of a dynamic system with x t the state of the system at time t and y t the observation at time t. The SS model can where, (a) the state vector is x t ∈ R nx with n x the state dimension, (b) the state transition density g(x t |x t−1 , θ) for t ≥ 2 and g(x 1 |x 0 , θ) := µ(x 1 |θ) describes the evolution of the states, (c) the observation is y t ∈ Y ny with n y the dimension of the observations, (d) the observation density conditional on the state is f (y t |x t , θ), and (e) the model includes a vector of unknown parameters θ ∈ Θ. The observations y 1:T = {y 1 , y 2 , ..., y T }, corresponding to the states x 1:T = {x 1 , x 2 , ..., x T }, are the time series data that can be stock returns, electro-magnetic signals, camera images, etc.
State filtering and likelihood estimation can be performed in a sequential Bayesian simulation framework by approximating p(x 1 |y 1 , θ) and p(y 1 |θ) at time t = 1, p(x 2 |y 1:2 , θ) and p(y 1:2 |θ) at time t = 2, and so on. This is known as the PF in the literature. In the standard PF, it is necessary that we can generate from the state transition density g(x t |x t−1 , θ) and can evaluate the observation density f (y t |x t , θ). Auxiliary particle filters often require that we can evaluate the state transition density g(x t |x t−1 , θ) as well as the observation density.
For the rest of the paper, we do not show dependence on θ, unless it is required.

Disturbance SS model
The PF for the conventional SS model described in section 2 assumes that it is possible to either generate from or evaluate the state transition density g(x t |x t−1 ). In some applications this is not possible or it is computationally expensive to compute (Murray et al., 2013;Hall et al., 2014). Suppose that, for t ≥ 2, we can write is a deterministic function and u t ∈ R nu is a vector of latent variables with density p U (u t ), such that k(u t , x t−1 )|x t−1 ∼ g(x t |x t−1 ) when u t ∼ p U (u t ). For t = 1, k(u 1 , x 0 ) := κ(u 1 ) with a deterministic function κ(·) such that κ(u t ) ∼ µ(x 1 ) when u 1 ∼ p U (u 1 ). Typically, u t is a set of independent uniform or standard normal random variables. We follow Murray et al. (2013) and refer to the u t as disturbances. Murray et al. (2013) propose reformulating the SS model in terms of the disturbances variables u 1:T where x t = x t (u 1:t ) = k(u t , x t−1 ). Murray et al. (2013) call this the disturbance SS model. Figure 2 gives its graphical dependence. It is possible to use a PF for the disturbance SS model (2) as the state transition density p(u t |u t−1 ) = p U (u t ) 6 is tractable. We note that, unlike the conventional SS model (1), the conditional distribution of y t in (2) depends on all the state disturbances u 1:t up to time t. State filtering and likelihood estimation proceeds by estimating p(u 1 |y 1 ) and p(y 1 ) at time t = 1, p(u 1:2 |y 1:2 ) and p(y 1:2 ) at time t = 2, etc. The likelihood in the conventional Figure 2: Pictorial representation of the disturbance SS model. (1) is

SS model
and the likelihood in the disturbance SS model (2) is Proposition 1 shows that (3) and (4) are equal. Its proof is obvious and omitted.
Proposition 1. Suppose that for each value of x t−1 , the mapping u t → x t = k(u t , x t−1 ) from R nu into R nx is one-to-one. Then (3) and (4) are equal.
Therefore we can carry out Bayesian inference for the model parameters θ using the PM approach, where the likelihood (4) is estimated unbiasedly based on a PF operating on the disturbance SS model (2). The next section presents the PF methodology for the disturbance SS model. 7 3.1 Particle filter for the disturbance SS model Sequential particle filtering for the disturbance SS model involves approximating p(u 1:t |y 1:t ) and the likelihood L t =´p(u 1:t , y 1:t )du 1:t , based on a set of weighted .. The conventional state can be recovered at every time step by x i t = k(u i t , x i t−1 ). The PF scheme is basically a sequential importance sampling and resampling procedure with the proposal density at step t having the following structure q t (u 1:t ) = q t−1 (u 1:t−1 )q t (u t |u 1:t−1 ).
Note that q t (u t |u 1:t−1 ) might depend on the data y 1:t as in the fully adapted PF.
The unnomarlized weights w t (u 1:t ) = p(u 1:t , y 1:t )/q t (u 1:t ) are decomposed as It is easy to show that L t L t−1 =´α t (u 1:t )w t−1 (u 1:t−1 )q t (u 1:t )du 1:t w t−1 (u 1:t−1 )q t (u 1:t )du 1:t , which can be estimated by with W i t−1 the normalized weights from time t − 1. Therefore the product of the estimates in (5) up to time t forms an estimate L t of the likelihood L t . By definition, L 0 = 1. It is well-known that L t is an unbiased estimate of L t (Del Moral, 2004).
Algorithm 1 gives the pseudo-code of the disturbance PF.
1:t and reset weights W i t = 1/N }

CPM for the disturbance SS model
The PF provides an unbiased estimate of the likelihood and this property facilitates the development of estimation techniques using the PM as in Andrieu et al. (2010), Murray et al. (2013) and Hall et al. (2014). This section proposes a CPM algorithm for Bayesian inference in the disturbance SS model.
The unbiased non-negative PF likelihood estimator can be written as L(θ, U), a function of the model parameters and the collection U of all the standard normal random numbers (or equivalently uniform random numbers) used to obtain the likelihood estimate. The collection U consists of, in the conventional PF, the N n x −dimensional random vectors used to generate the new states x t at each time step t = 1, ..., T , and in the disturbance PF, the N n u −dimensional random vectors used to generate u t at each time step t = 1, ..., T . Moreover resampling involves additional random numbers. This paper employs systematic resampling which requires one random number per resampling step. Denote by p U (U) the density of U. The CPM method is derived as follows. Define the joint pseudo target of θ and U as The marginal π(θ) of π(θ, U) with respect to θ is the same as the posterior π(θ) = p(y 1:T |θ)p Θ (θ)/p(y 1:T ), because´ L(θ, U) p U (U)dU = p(y 1:T |θ). That is, we can run an MCMC on the expanded space of (θ, U), and obtain iterates from π(θ). Let θ c and U c be the current values of the parameters and the random numbers respectively.
To iterate through the MCMC chain, we first propose their corresponding new values according to and then accept this proposal with the probability α = min 1, This standard PM method has been applied to both the conventional and the disturbance SS models (Andrieu et al., 2010;Pitt et al., 2012;Murray et al., 2013).
An important version of the PM algorithm is the CPM algorithm of Dahlin et al.
(2015a) and Deligiannidis et al. (2016) where the likelihood estimates appearing in the MH ratio are correlated by correlating the current and proposed standard normal random numbers U as where ρ is the non-negative correlation between the random numbers and ξ is a standard normal vector of the same size as U. The standard PM is the special case of CPM when ρ = 0. The CPM algorithm accelerates the estimation process significantly. The strong correlation between logs of the likelihood estimates reduces the variability in the likelihood ratio of the proposed and current parameter values resulting in faster convergence.
In this paper, we propose using CPM for the disturbance SS model. It is important that the correlation in the random numbers is preserved in the estimated likelihoods. However the resampling step of the PF might impede this preservation because of its particle replacement property (Deligiannidis et al., 2016). A small change in the random numbers used in the resampling steps might lead to a big change in the particle paths, and thus the correlation in the logs of the likelihood estimates might not be well preserved. This impediment can be facilitated in two steps.
Fixed resampling frequency: Resampling is usually performed whenever the effective sample size N eff falls below a threshold. In our CPM the particles are resampled once every R f time steps with 1 ≤ R f ≤ T . Having a deterministically fixed resampling frequency R f allows particles to be prearranged in a way that preserves the correlation.
Particle ordering: Prior to being resampled, the particles are ordered so that they are close to one another in some metric. Ordering can be easily accomplished in the univariate case by sorting the particles from smallest to the largest. For the multivariate case, such sorting is unavailable. One approach for ordering multidimensional particles is to use a Hilbert space filling curve method (Skilling, 2004), which is used in the PF by Gerber and Chopin (2015) (in a context not related to CPM), and is proposed for the CPM algorithm by Deligiannidis et al. (2016). The Hilbert curve method transforms multidimensional particles to a univariate space (the Hilbert space in this context) based on some metric so that particle locality is preserved. The resulting transformed univariate set of particles is then sorted to obtain the sorting indices, i.e., the Hilbert curve method provides a mechanism to transform multidimensional particles to a space on which traditional sorting can be applied. Our article proposes the following simpler and more resource-efficient multidimensional sorting scheme.

Multidimensional Euclidean sorting
..,N be the n x −dimensional particles at a given time step, The time subscript is removed for notational simplicity. Let d(x j , x i ) be the Euclidean distance between two multidimensional particles x i and x j . Algorithm 2 describes the procedure to generate the set of sorting indices for the particles. The first sorting index in the algorithm is the index of the particle having the least value along its first dimension. The rest are chosen in a way that minimizes the Euclidean distance between the recently selected particle and the set of all remaining particles.
This approach is employed in our article to sort multidimensional particles.
The two PF modifications -fixed resampling frequency and particle orderingensure that the variation in the particle states due to resampling is minimized. Hence the correlation in the likelihood is more likely to be preserved leading to good mixing of the MCMC chain with fewer particles.

CPM for the disturbance SS models
The CPM can be extended to the disturbance SS models by ordering the original states. However, we take a different approach by working in the disturbance space.
The justification is that, since the particles are often a smooth function of the disturbances, a similarity (or variation) in the disturbances leads to similarity (or variation) in the state particle. If the disturbances are close to each other then the corresponding state particles will also be close to each other. This implies that sorting the disturbances is sufficient to preserve the correlation in the PF likelihood values. The

BPM for SS models
The CPM method involves particle sorting in order to maintain proximity within the particles. This procedure can be computationally expensive and renders the CPM method infeasible for high dimensional models involving a large number of observations. Furthermore, the correlation in the log likelihood estimates is controlled by the correlation ρ between the U, but this relation is model-dependent and not known precisely. A large ρ does not necessarily lead to a high correlation in the log likelihood estimates (see section 6). These limitations can be overcome in the BPM approach proposed in this section.
The BPM method of Tran et al. (2016) is an alternative to the CPM, in which the set U is divided into G blocks, each of these is updated jointly with θ in each MCMC iteration. Under some assumptions, they show that the optimal number of particles should be selected such that the variance of the log-likelihood estimate is σ 2 opt = 2.16 2 /(1 − ρ 2 l ) with ρ l the correlation between the log of the likelihood estimates. We extend the BPM to the general class of SS models where the unbiased likelihood estimate is an average of G likelihood estimates obtained by G independent PFs.
Let L θ, U (i) be the estimated likelihood obtained from the ith PF, i = 1, ..., G.
We define the joint target density of θ and U = U (1) , . . . , U (G) as is the average of the G unbiased likelihood estimates and hence also unbiased. We then update the parameters jointly with a randomly-selected block U (K) in each MCMC iteration, with Pr (K = k) = 1/G for any k = 1, ..., G. Using this scheme, the acceptance probability is The PFs can be run in parallel on multiple processors. It is possible to show that the correlation between log L (θ p , U p ) and log L (θ c , U c ) is approximately ρ l = 1 − 1/G, so the more particle filters we run the higher the correlation. Unlike the CPM, the BPM method allows a more direct control of the correlation between the logs of the where Time is the elapsed time (in seconds) per iteration. We also compute the relative time normalized variance (RTNV) defined as where the benchmark is chosen to be the standard PM method for all the models studied in this paper.
In all the studies, the total number of MCMC iterations is 25000, with the first 5000 discarded as burnin iterations. The particles are resampled at every time step.
The adaptive random walk (ARW) algorithm of Roberts and Rosenthal (2009) is used, with the scale in the ARW proposal adapted using the Robbins-Monro method of Garthwaite et al. (2016) to obtain an acceptance rate of 0.23. The proposal using the gradient information from the likelihood is also used for the BPM method. See the Appendix for the details. The number of particles in the standard PM method is chosen such that Var(log L(θ)) ≈ 1. The number of particles in the CPM methods (with particle and disturbance sorting) is chosen such that Var(log L(θ)) ≈ 2.
where ρ l is the log likelihood correlation estimated over 50 independent replications of the particle filter with the parameters fixed at their true values. For the BPM method, we use G = 12 PFs, hence ρ l = 1 − 1/12 = 0.9167 and the number of particles is chosen such that Var(log L(θ)) ≈ 29.21.

Stochastic Volatility model
The The state transition and the observation equations of the corresponding SS model The transition equation can be written as u t ∼ N (0, 1). The original state dimensionality is n x = P and the disturbance dimensionality is n u = 1, hence when P > 1 the disturbance has a lower dimension than the state.
We study the AR(4) SV model to demonstrate the efficiency gain achieved by the proposed CPM method operating in the univariate disturbance space over the CPM method operating in the four-dimensional conventional SS. The simulation was carried out for two cases with T = 1000 and T = 3000. The CPM methods (with particle sorting) was omitted in the T = 3000 case because of the computational complexity involved in sorting the multidimensional particles. We note that the estimated correlation of the logs of the likelihood in the CPM method decreases with the dimension of the state: a random number correlation ρ = 0.9999 induced only ρ l ≈ 0.7 correlation within the log likelihood estimates in the AR(4) SV model, while for a univariate SV model we obtain ρ l > 0.99.
Our results show that the proposed CPM operating in the disturbance space is able to target almost the same log likelihood variance as targeted by the standard CPM method operating on the original state space. That is, the proposed CPM method uses as many particles as required by the standard CPM method and as few as 15% of that required by the standard PM method. This validates the key claim of this paper that sorting the disturbances is sufficient to maintain order within the particles in the way that it preserves the log likelihood correlation as much as that preserved by directly sorting the multidimensional particles. Moreover the MCMC convergence obtained by CPM method using the proposed multidimensional Euclidean sorting is commensurate with that of the Hilbert curve method. Hence the proposed Euclidean sorting scheme is valid and applicable to the CPM scheme. It can also be observed from the table that the proposed CPM method for the disturbance SS models is ≈ 2 times and ≈ 9 times better than the standard PM in terms of RTNV for T = 1000 and T = 3000 respectively. A tremendous gain over the conventional CPM method using the Hilbert curve and Euclidean sorting can also be observed: ≈ 12 times better than CPM with Euclidean sorting and ≈ 30 better than CPM with Hilbert sorting at T = 1000. This improvement will increase with the target state dimensionality and the number of observations. It is noteworthy that although the conventional CPM method operating in the target space facilitates the use of fewer particles for parameter estimation, the method does not sustain its gain in terms of RTNV, at small values of T, due to the inclusion of the resource hungry multidimensional sorting procedure. This advocates the need for designing resource efficient multidimensional sorting methods or improved filtering schemes, e.g., the efficient importance sampling Jean-Francois and Zhang (2007)

Nonstationary growth model
The univariate nonlinear nonstationary growth model is described by the SS model where x 1 ∼ N (0, 5), a t ∼ N (0, τ 2 ), and e t ∼ N (0, σ 2 ). The parameters for this model are θ = (τ 2 , σ 2 ). This is a popular example used in the literature to assess the performance of the particle filter for nonlinear models. The number of observations is set to T = 200. The true parameter values are τ 2 true = 10 and σ 2 true = 1. The priors for the parameters log τ 2 ∼ N (0, 50) and log σ 2 ∼ N (0, 50). Table 2 summarizes the estimation results. Figure 5 shows the corresponding trace plots and the ACFs We observe that the four PM samplers converge adequately to the same invariant distribution.
Note that in the CPM method we were unable to target the optimal variance. Our tests revealed that, with the number of particles N chosen such that Var(log L(θ)) ≈ 2.16 2 / (1 − ρ 2 l ) with ρ l = 0.9868, the MCMC chain is sticky and the IACT value is large. This behavior is observed consistently in nonstationary models with multivariate states, e.g., the spline model, the bearings only tracking model, etc. Hence in our study, we increased the number of particles so that the MCMC chain mixes well.
In contrast, our tests revealed that the proposed BPM method applies to a general class of SS models and induces a precisely known and controllable (via the number of parallel cores) correlation within the log likelihood estimates.
The table shows that in terms of the TNV, the BPM method is ≈ 2.5 times better than CPM and the standard PM methods. Note that including gradient information in the proposal improves the efficiency of the BPM sampler in terms of the IACT and the method is ≈ 1.5 times better than the ARW BPM method. It is possible that increasing the number of particles will make the gradient more informative (i.e.,  reduce its variance) and will hence improve the parameter proposal.

Cubic spline model
A simple bivariate cubic spline model is defined by the SS formulation where the state transition matrix is F (δ t−1 ) =    N (0, σ 2 ). The rate at which observations are received is chosen to be δ t = 1/T .
The number of observations is set to T = 500. The true parameter values are τ 2 true = 4 and σ 2 true = 0.25. The priors for the parameters are τ 2 ∼ IG(1, 1) and 1). The BPM method with gradient information is not included in this study because the gradient involves the inversion of the matrix U (δ t−1 ), which has a large value, thus making the gradient computation numerically unstable. Table 3 summarizes the estimation results. Figure 6 shows the corresponding trace plots and the ACFs. The figure shows that the PM samplers converge adequately to the same invariant distribution and that the estimated mean and standard deviation values are similar. The results also show that the BPM method also exhibits improved performance in terms of IACT and is ≈ 7 times better than the standard PM method.

Lotka Volterra model
The Lotka Volterra model is a reaction network that models the interaction between species (Golightly and Wilkinson, 2011). Here we study a simple bivariate nonstationary stochastic Lotka Volterra model. Define the state vector by x = (x 1 , x 2 ) , where x 1 denotes the prey and x 2 the predators. The predator prey system is comprised of three reactions  representing the prey reproduction, the predator-prey interaction and the predator death. The rate constant vector for the three reactions is defined as c = (c 1 , c 2 , c 3 ) and the associated hazard function is h(x, c) = (c 1 x 1,t , c 2 x 1,t x 2,t , c 3 x 2,t ) .
The stoichiometry matrix used to update the state if a reaction occurs is given We set x 1 = (100, 100) and the (intractable) forward simulations are conducted using the Gillespie algorithm Wilkinson (2011). The observation equation is where the noise variance e t ∼ N (0, σ 2 ). The parameters for this model are θ = {c, σ 2 }.
The number of observations is T = 50. The true parameter values are c true = (0.5, 0.0025, 0.3) and σ 2 true = 0.5. The priors for the parameters are c i=1:3 ∼ U(0, 1) and σ 2 ∼ IG(1, 1). Table 4 summarizes the estimation results. Figure 7 plots the corresponding trace plots and the ACFs, and shows that the two PM samplers converge adequately to the same invariant distribution. The results suggest that the BPM method exhibits improved performance in terms of IACT and is ≈ 40 times better than the standard PM method.

Application of the stochastic volatility model to financial returns data
This section studies the performance of the proposed PM methods in financial data.
We consider two data sets corresponding to, (a) the daily returns for the Australian The first data set is the All Ords index and contains T = 3564 observations. The second data set corresponds to the S&P500 index and contains T = 7832 observations.
The data is retrieved from Yahoo finance. We fit the AR(4) SV model described in section 6.1 to these two data sets. The parameters for the models are θ = (φ, τ 2 ).
The priors for the parameters and the ARW proposal are the same as described in section 6.1. Table 5 reports the estimation results corresponding to the All Ords data set, for the standard PM, disturbance CPM and the BPM methods. The BPM method was implemented on G = 12 independent cores. The result corresponding to the CPM with particle sorting is omitted due to the computational infeasibility of the method. Figure 8 shows the corresponding trace plots and the ACFs, and shows that the PM samplers converge adequately to the same invariant distribution. The disturbance CPM method uses 4.3 times fewer particles and is 26 times superior to the standard PM method. The clear winner is the BPM method. The BPM (with an ARW proposal) method uses 26 times fewer particles than the standard PM method, and is ≈ 838 times and ≈ 33 times superior to the standard PM and the disturbance CPM methods respectively. Note that the number of particles used in the BPM with the proposal using the derivative information is slightly increased to make the gradient more informative. This leads to reduced IACT values. It can be observed that the BPM method (with the proposal using the Hessian information) is ≈ 1540 times and ≈ 57 times superior to the standard PM and the disturbance CPM methods respectively.  T = 7832, for the disturbance based CPM and the BPM methods. The standard PM requires 3000 particles to satisfy Var(log L(θ)) ≈ 1. Consequently its TNV is large.
The results corresponding to the standard PM and the CPM with particle sorting are excluded due to the computational infeasibility of the methods. Figure 9 reports the corresponding trace plots and the ACFs and shows that the PM samplers converge adequately to the same invariant distribution. The BPM method is implemented on G = 4 independent cores. The BPM methods are ≈ 100 and ≈ 280 times superior to the disturbance CPM method using only half the number of particles used in the CPM. The BPM method using the gradient information are ≈ 3 times better than the ARW BPM method.

Conclusion
The contribution of this paper is two-fold, (a) a CPM method for parameter estimation in disturbance SS models is proposed, and (b) a BPM method for general time series models is proposed. The main insight in the CPM method for disturbance SS models is that the disturbance SS possesses the same properties as that the conventional SS as required to preserve the correlation within the likelihood estimates.
This insight is the basis for the proposed CPM method for disturbance SS models.
The key innovation of this proposal is that the CPM is performed in the disturbance space rather than the SS by virtue of preserving the likelihood correlation by maintaining similarity in the disturbance states as opposed to the traditional approach that maintains similarity in the target states. The main insight in the BPM method is that implementing multiple independent particle filters, with the likelihood estimated by averaging the likelihood estimates obtained from separate particle filters, facilitates the reduction in the variance of the log of the estimated likelihood and hence in the use of fewer particles. This insight is the basis for the design of a flexible means to block the random numbers in order to obtain a desired correlation for the log of likelihood estimates. A Proposal using information from the derivatives This section presents the construction of the proposal density in MCMC that makes use of the derivatives of the log likelihood. Poyiadjis et al. (2011) were the first to show how the particle filter methods can be used to estimate the derivatives of the log likelihood for state space models. Their methods might suffer from a computational cost that is quadratic in the number of particles, Nemeth et al. (2016) proposed an alternative method whose computational cost is linear in the number of particles.
They use a combination of kernel density estimation and Rao-Blackwellisation to reduce the Monte Carlo error of the estimates. For non-linear and non-Gaussian state space models it is impossible to obtain the exact derivatives, but they can be approximated using the particle approximation of p (x 1:T |y 1:T , θ).
Setting λ = 1 gives the algorithm in Poyiadjis et al. (2011). Nemeth et al. (2016) show that the bias and variance of both score estimate and observed information matrix vary according to λ. Reducing the value of λ increases the bias, but it reduces the Monte Carlo variance of estimates. It is also shown that setting λ ≈ 0.95 will produce an estimate for the score and observed information matrix with linearly increasing variance and minimal bias. Therefore we use λ = 0.95 in all our applications. The parameter proposal adopted in this paper is similar to the one used in Dahlin et al.
where G (θ; U) = 1 Algorithm 3 Algorithm to estimate Gradient and Hessian Matrix • At iteration t = 1, ..., T -Run the particle filter to obtain x t is the weight of particle i at time t. k i is the ancestor index of particle i at time t − 1.
-Normalise the weights W • The observed negative Hessian matrix can be estimated as  Table 1. The columns from left to right correspond to, (a) the standard PM, (b) the CPM with sorting the multidimensional particles using the Hilbert curves method, (c) the CPM with sorting the multidimensional particles using the proposed Euclidean sorting method, (d) the CPM with sorting the univariate disturbances, (e) the proposed BPM method, and (f) the BPM method with the proposal using the derivative information.  Figure 8: The estimation result for the Australian All Ords financial returns data applied to the AR(4) SV model corresponding to the result reported in Table 5. The columns from left to right correspond to, (a) the standard PM, (b) the CPM with sorting the univariate disturbances, (c) the proposed BPM method with RW, and (d) the proposed BPM method with the proposal using the derivative information.  Figure 9: The estimation result for the S&P 500 daily financial returns data applied to the AR(4) SV model corresponding to the result reported in Table 5. The columns from left to right correspond to, (a) the CPM with sorting the univariate disturbances, (b) the BPM with RW, and (c) the BPM with the proposal using the derivative information.