Particle rejuvenation of Rao-Blackwellized Sequential Monte Carlo smoothers for Conditionally Linear and Gaussian models

# Particle rejuvenation of Rao-Blackwellized Sequential Monte Carlo smoothers for Conditionally Linear and Gaussian models

Ngoc Minh Nguyen111LTCI, CNRS and Télécom ParisTech, 46 rue Barrault 75634 Paris Cedex 13, France.    Sylvain Le Corff222Laboratoire de Mathématiques d’Orsay, Univ. Paris-Sud, CNRS, Université Paris-Saclay, 91405 Orsay, France.    Eric Moulines333Centre de Mathématiques Appliquées, UMR 7641, Ecole Polytechnique, France.
###### Abstract

This paper focuses on Sequential Monte Carlo approximations of smoothing distributions in conditionally linear and Gaussian state spaces. To reduce Monte Carlo variance of smoothers, it is typical in these models to use Rao-Blackwellization: particle approximation is used to sample sequences of hidden regimes while the Gaussian states are explicitly integrated conditional on the sequence of regimes and observations, using variants of the Kalman filter / smoother. The first successful attempt to use Rao-Blackwellization for smoothing extends the Bryson-Frazier smoother for Gaussian linear state space models using the generalized two-filter formula together with Kalman filters / smoothers. More recently, a forward backward decomposition of smoothing distributions mimicking the Rauch-Tung-Striebel smoother for the regimes combined with backward Kalman updates has been introduced. This paper investigates the benefit of introducing additional rejuvenation steps in all these algorithms to sample at each time instant new regimes conditional on the forward and backward particles. This defines particle based approximations of the smoothing distributions whose support is not restricted to the set of particles sampled in the forward or backward filter. These procedures are applied to commodity markets which are described using a two factor model based on the spot price and a convenience yield for crude oil data.

## 1 Introduction

State space models are bivariate stochastic processes where the state sequence is a Markov chain which is only partially observed through the sequence . Conditionally on the state sequence the observations are independent and for all the conditional distribution of given depends on only. These models are used in a large variety of disciplines such as financial econometrics, biology, signal processing, see [DM13] and the references therein. In general state space models, bayesian filtering and smoothing problems, i.e. the computation of the posterior distributions of a sequence of states for given observations , are challenging tasks. Filtering refers to the estimation of the distributions of the hidden state given the observations up to time , while fixed-interval smoothing stands for the estimation of the distribution of sequence of states given observations with .

When the state and observation models are linear and Gaussian, filtering can be solved explicitly using the Kalman filter [Kal60]. Exact solutions of the fixed-horizon smoothing problem can be obtained using either the Rauch-Tung-Striebel smoother [RST65] or the Bryson-Frazier two-filter smoother [BF63]. This paper focuses on Conditionally Linear and Gaussian Models (CLGM) given for by:

 Zi=dai+TaiZi−1+Haiεi, (1)

where:

1. is a sequence of independent and identically distributed (i.i.d.) -dimensional Gaussian vectors with zero mean and identity covariance.

2. is a homogeneous Markov chain taking values in a finite space , called regimes, with initial distribution and transition matrix .

3. are positive-definite matrices, -dimensional vectors and positive-definite matrices.

4. is a -dimensional Gaussian random variable with mean and variance independent of .

Let be the number of observations. At each time step , the observation is given by:

 Yi=cai+BaiZi+Gaiηi, (2)

where:

1. is a i.i.d. sequence of -dimensional Gaussian vectors, independent of and .

2. are positive-definite matrices, -dimensional vectors and matrices.

CLGM play an important role in many applications; see [Sar13] and the references therein for an up-to-date account. A crucial feature of these models is that, conditional on the regime sequence , both the state equation and the observation equation are linear and Gaussian, which implies that conditional on the sequence of regimes and on the observations, the filtering and the smoothing distributions of the continuous states can be computed explicitly.

To exploit this specific structure, it has been suggested in the pioneering works of [CL00, DGA00] to solve the filtering problem by combining Sequential Monte Carlo (SMC) methods to sample the regimes with the Kalman filter to compute the conditional distribution of the states sequence conditional on the regimes and on the observations. This is a specific instance of Rao-Blackwellized Monte Carlo filters, often referred to as the Mixture Kalman Filter. Improvements of these early filtering techniques have been introduced in [DGK01, SGN05].

The use of Rao-Blackwellization to solve the smoothing problem has been proved to be more challenging and has received satisfactory solutions only recently. The first forward-backward smoother proposed in the literature [FGDW02] was not fully Rao-Blackwellized as it required to sample the hidden linear states in the backward pass. An alternative approach, based on the so-called structural approximation of the model suggested in an early paper by [Kim94], was proposed by [Bar06] to avoid to sample a continuous state in the backward pass. This approximation is rather ad-hoc and the resulting smoother is not consistent when the number of particles goes to infinity. The inaccuracy introduced by the approximation might be difficult to control.

The first fully Rao-Blackwellized SMC smoother which should lead to consistent approximations when the number of particles grows to infinity was proposed by [BDM10] and extends the Bryson-Frazier smoother for Gaussian linear state space models using the generalized two-filter formula with Rao-Blackwellization steps for the forward and the backward filters. This two-filter approach combines a forward filter with a backward information filter which are approximated numerically using SMC for the regime sequence and Kalman filtering techniques for the hidden linear states.

More recently, [SBG12, LBGS13, LBS16] introduced a Rao-Blackwellized smoother based on the forward-backward decomposition of the FFBS algorithm with Rao-Blackwellization steps both in the forward and backward time directions. The update of the smoothing distribution of the regime given the observations shares some striking similarities with the Rauch-Tung-Striebel smoothing procedure, which is at the heart of the FFBS procedure. The Rao-Blackwellization requires to update backward in time the smoothing distribution of the states given the regimes and the observations, which is achieved by using an à la Kalman backward update.

In this paper, we propose to improve the performance of the algorithms introduced in [BDM10] and in [SBG12, LBGS13, LBS16] by using additional Rao-Blackwellization steps which allows to sample new particles in the backward pass. This approach may be seen as an extension of the ideas of [FC03] for Rao-Blackwellized smoothers. In [BDM10], for all , the sampled forward and backward sequences are merged to approximate the posterior distribution of . This provides an approximation whose support is restricted to the particles produced at time by the backward particle filter. As noted in [FWT10, Secion 2.6], these two-filter smoothers are prone to suffer from degeneracy issues when the algorithm associates forward particles at time with backward particles at time . We propose to approximate the marginal smoothing distribution of by merging the sampled forward and backward trajectories at times and and integrating out all possible paths between time and time and between time and time instead of sampling random variables. Similarly, in the backward pass of [SBG12, LBGS13, LBS16], a regime is sampled at time using the particles produced by the forward filter at time . In this case, particle rejuvenation may be introduced by using the forward weighted samples at time and extending these trajectories at time with a Kalman filter for all possible values of the regime. Then, is sampled in using an appropriately adapted weight.

The paper is organized as follows. The algorithms introduced in [BDM10] and in [SBG12, LBGS13, LBS16] as long as the proposed rejuvenation associated with each method are presented in Section 2. The performance of all these methods is illustrated in Section 3 with simulated data. In Section 4, an application to commodity markets is presented; the performance of our procedure is illustrated with crude oil data. A detailed derivation of the algorithms is provided in the Appendix.

## 2 Rao-Blackwellized smoothing algorithms

This section details the Sequential Monte Carlo algorithms which can be used to approximate the conditional distribution of the states or the marginal distributions of given the observations . For all matrix let be the determinant of . If is a positive-definite matrix, for all define

 ∥z∥2A:=z′A−1z,

where for any vector or matrix , denotes the transpose matrix of . Let be the probability density of the conditional distribution of given and be the probability density of the conditional distribution of given :

 m(ai,zi−1;zi) :=∣∣2π¯¯¯¯¯Hai∣∣−1/2exp{−12∥∥zi−dai−Taizi−1∥∥2¯¯¯¯Hai}, (3) g(ai,zi;yi) :=∣∣2π¯¯¯¯Gai∣∣−1/2exp{−12∥∥yi−cai−Baizi∥∥2¯¯¯¯Gai}, (4)

where

 ¯¯¯¯Gj:=GjG′j,¯¯¯¯¯Hj:=HjH′j.

All the algorithms considered in this paper are based on forward-backward or two-filter decompositions of the smoothing distributions and share the same forward filter presented in Section 2.1.

### 2.1 Forward filter

The SMC approximation of may be obtained using a standard Rao-Blackwellized algorithm. The procedure produces a sequence of trajectories associated with normalized importance weights () used to define the following approximation of :

 pN(a1:i,zi|y1:i)=N∑k=1ωkip(zi|ak1:i,y1:i)δak1:i(a1:i), (5)

where is the Dirac delta function. In this equation, the conditional distribution of the hidden state given the observations and a trajectory is a Gaussian distribution whose mean and variance may be obtained by using the Kalman filter update.

Initialization
At time , write, for all ,

 μj1|0=cj+Bjμ1andPj1|0=BjΣ1B′j+¯¯¯¯Gj.

are sampled independently in with probabilities proportional to

 πjp(a1=j|y1)∝πj|Pj1|0|−1/2exp{−(y1−μj1|0)′(Pj1|0)−1(y1−μj1|0)/2}.

Then, and are computed using a Kalman filter:

 Kk1 =Σ1B′ak1(Bak1Σ1B′ak1+¯¯¯¯Gak1)−1, μk1 =μ1+Kk1(Y1−cak1−Bak1μ1), Pk1 =(Im−Kk1Bak1)Σ1,

where for all positive integer , is the identity matrix. Each particle particle is associated with the importance weight .

Iterations
Several procedures may be used to extend the trajectories at time . For all sampled trajectories and all , [CL00] used the incremental weights:

 γj,ki=p(yi|ai=j,ak1:i−1,y1:i−1)Q(aki−1,j).

The conditional distribution of given , and is a Gaussian distribution with mean and variance where

 μki|i−1(ai) =dai+Taiμki−1, Pki|i−1(ai) =TaiPki−1T′ai+¯¯¯¯¯Hai.

Therefore,

 γj,ki∝Q(aki−1,j)|BjPj,ki|i−1B′j+¯¯¯¯Gj|−1/2exp{−12∥∥yi−cj−Bjμj,ki|i−1∥∥2BjPj,ki|i−1B′j+¯¯¯¯Gj},

where:

 μj,ki|i−1 =μki|i−1(j)=dj+Tjμki−1, (6) Pj,ki|i−1 =Pki|i−1(j)=TjPki−1T′j+¯¯¯¯¯Hj. (7)

In [CL00], for all , an ancestral path is chosen with probabilities proportional to . Then, the new regime is sampled in with probabilities proportional to . A drawback of this method is that only ancestral paths that have been selected using the importance weights are extended at time . Following [BGM08], this may be improved by considering all the offsprings of all ancestral trajectories . Each ancestral path has offsprings at time , it is thus necessary to choose a given number of trajectories at time (for instance ) among the possible paths. To obtain the weight associated with each offspring write the following approximation of based on the weighted samples at time :

 pN(a1:i|y1:i) ∝N∑k=1ωki−1Q(aki−1,ai)p(yi|ak1:i−1,ai,y1:i−1)δak1:i−1(a1:i−1), ∝N∑k=1J∑j=1ωki−1γj,kiδ(ak1:i−1,j)(a1:i).

Therefore, each ancestral trajectory of the form , , , is associated with the normalized weight . Several random selection schemes have been proposed to discard some of the possible offsprings to maintain an average number of particles at each time step. Following [BGM08], we might choose between the Kullback-Leibler Optimal Selection (KL-OS) or the Chi-Squared Optimal Selection (CS-OS) to associate a new weight to each of the trajectories. If the new weight is 0, then the corresponding particle can be removed.

KL-OS: is chosen as the solution of :

 N∑k=1J∑j=1min(~ωj,ki/λ,1)=N.

For all and , if then the new weight is and if :

 ~Ωj,ki={λwith probability ~ωj,ki/λ,0with probability 1−~ωj,ki/λ.

CS-OS: is chosen as the solution of :

 N∑k=1J∑j=1min(√~ωj,ki/λ,1)=N.

For all and , if then the new weight is and if :

 ~Ωj,ki=⎧⎪ ⎪⎨⎪ ⎪⎩√~ωj,kiλwith probability √~ωj,ki/λ,0with probability 1−√~ωj,ki/λ.

Then, in both cases, all particles such that are discarded and for all the other trajectories defined as an ancestral path extended by , the new corresponding weight in (5) is given by the normalized weight . In the numerical sections of this paper, the Kullback-Leibler Optimal Selection (KL-OS) scheme is used.

### 2.2 FFBS based algorithms

#### 2.2.1 FFBS algorithms of [Sbg12, Lbgs13, Lbs+16]

[SBG12, LBGS13, LBS16] proposed a Rao-Blackwellized procedure to sample the regime backward in time following the same steps as in the Forward Filtering Backward Smoothing algorithm [HK98, DGA00]. The algorithm relies on the decomposition given, for all , by:

 p(a1:n|y1:n)=p(a1:i|ai+1:n,y1:n)p(ai+1:n|y1:n).

This decomposition is similar to the Rauch-Tung-Striebel decomposition of the filtering distribution. The first factor on the right hand side of the previous equation is nevertheless more difficult to handle because it itself relies on all the observations. As noted by [SBG12], this term can be computed recursively by considering the following decomposition:

 p(a1:i|ai+1:n,y1:n)∝p(yi+1:n,ai+1:n|a1:i,y1:i)p(a1:i|y1:i). (8)

The second factor in the last equation may be approximated using the ancestral trajectories and the associated importance weights produced by the forward filter. Therefore, may be approximated by:

 pN(a1:i|ai+1:n,y1:n)=N∑k=1~ωki|nδak1:i(a1:i)with~ωki|n∝ωkip(yi+1:n,ai+1:n|ak1:i,y1:i).

Then, a trajectory approximatively distributed according to may be sampled following these steps:

1. Set with probabilities proportional to .

2. For all , set with probabilities proportional to .

This algorithm requires to compute the quantity . This predictive quantity is available analytically using Kalman filtering techniques. However, this has to be done for each trajectory , which leads to an algorithm with a prohibitive computational complexity. [LBS16] proposed a procedure computationally less intensive by conditioning with respect to and then marginalizing with respect to this variable:

 p(yi+1:n,ai+1:n|ak1:i,y1:i)=∫p(yi+1:n,ai+1:n|zi,aki)p(zi|ak1:i,y1:i)dzi. (9)

This is similar to the two-filter decomposition of the smoothing distribution, see Section 2.3. By [LBS16],

 p(yi+1:n,ai+1:n|zi,ai)∝Q(ai,ai+1)exp{−(z′iΩi(ai+1:n)zi−2λ′i(ai+1:n)zi)/2},

where the proportionality is with respect to and

 p(yi:n,ai+1:n|zi,ai)∝exp{−(z′iˆΩi(ai:n)zi−2ˆλ′i(ai:n)zi)/2},

where the proportionality is with respect to . These quantities may be computed recursively backward in time with:

 ˆΩn(an) =B′an¯¯¯¯G−1anBan, ˆλn(an) =B′an¯¯¯¯G−1an(yn−can).

Then, for , define and and write

 Ωi(ai+1:n) =T′ai+1(Im−ˆΩi+1(ai+1:n)Hai+1M−1i+1H′ai+1)ˆΩi+1(ai+1:n)Tai+1, λi(ai+1:n) =T′ai+1(Im−ˆΩi+1(ai+1:n)Hai+1M−1i+1H′ai+1)mi+1.

As ,

 ˆΩi(ai:n) =Ωi(ai+1:n)+B′ai¯¯¯¯G−1aiBai, ˆλi(ai:n) =λi(ai+1:n)+B′ai¯¯¯¯G−1ai(yi−cai).

Then, by (9),

 p(yi+1:n,ai+1:n|ak1:i,y1:i)∝Q(aki,ai+1)|Λki(ai+1:n)|−1/2exp{−ηki(ai+1:n)/2}, (10)

where the proportionality is with respect to and

 Λki(ai+1:n) =(Γki)′Ωi(ai+1:n)Γki+Im, ηki(ai+1:n) =∥μki∥2Ω−1i(ai+1:n)−2λ′i(ai+1:n)μki−∥(Γki)′(λi(ai+1:n)−Ωi(ai+1:n)μki)∥2Λi(ai+1:n),

where . Therefore,

 ~ωi|n∝ωkiQ(aki,ai+1)|Λki(ai+1:n)|−1/2exp{−ηki(ai+1:n)/2}.

If are independent copies of , the SMC approximation of [LBS16] of the joint smoothing distribution of the regime is:

 pLbscg~N(a1:n|Y1:n)=1~N~N∑k=1δ~ak1:n(a1:n).

#### 2.2.2 Particle rejuvenation of FFBS algorithms

The crucial step of the FFBS algorithm is the decomposition (8) which allows to extend a backward trajectory by choosing a particle in the set produced by the forward filter (and discarding the states ). An improved version of this FFBS algorithm which is not constrained to sample states in the support may be defined for all by writing:

 p(a1:i|ai+1:n,y1:n) ∝p(yi+1:n,ai+1:n|a1:i,y1:i)p(a1:i|y1:i), ∝p(yi+1:n,ai+1:n|a1:i,y1:i)∫p(a1:i−1,zi−1|y1:i−1)Q(ai−1,ai) m(ai,zi−1;zi)g(ai,zi;yi)dzi−1:i.

Replacing in the integral by the particle approximation obtained during the forward pass and using Kalman filtering techniques for each trajectory and each yields:

 ∫pN(a1:i−1,zi−1|y1:i−1)Q(ai−1,ai)m(ai,zi−1;zi)g(ai,zi;yi)dzi−1:i∝N∑k=1ωki|i−1(ai)δak1:i−1(a1:i−1),

where

 ωki|i−1(ai)=ωki−1Q(aki−1,ai)|Σki|i−1(ai)|−1/2exp{−12∥yi−yki|i−1(ai)∥Σki|i−1(ai)},
 yki|i−1(ai)=cai+Bai(dai+Taiμki−1)andΣki|i−1(ai)=Bai(TaiPki−1T′ai+¯¯¯¯¯Hai)B′ai+¯¯¯¯Gai.

On the other hand, for all , is computed as in (10) with all possible values and not only the regime of the filtering pass . This means that a Kalman filter must be used for each trajectory which may be extended by . Denote by and the mean and covariance matrix of the law of given obtained as in (6) and (7). Then,

 p(yi+1:n,ai+1:n|ak1:i−1,ai,y1:i)=Q(ai,ai+1)|Λki|i−1(ai:n)|−1/2exp{−ηki|i−1(ai:n)/2}, (11)

where the proportionality is with respect to and

 Λki|i−1(ai:n) =(Γki|i−1(ai))′Ωi(ai+1:n)Γki|i−1(ai)+Im, ηki|i−1(ai:n) =∥μki|i−1(ai)∥2Ω−1i(ai+1:n)−2λ′i(ai+1:n)μki|i−1(ai) −∥(Γki|i−1(ai))′(λi(ai+1:n)−Ωi(ai+1:n)μki|i−1(ai))∥2Λi(ai+1:n),

where is defined as . The distribution is then approximated by :

 pN(a1:i|ai+1:n,y1:n)∝N∑k=1ωki|i−1(ai)Q(ai,ai+1)|Λki|i−1(ai:n)|−1/2exp{−ηki|i−1(ai:n)/2}δak1:i−1(a1:i−1). (12)

By integrating over all possible paths , is sampled in . This particle rejuvenation step allows to explore states which are not in the support of the particles produced by the forward filter and improves the accuracy and the variance of the original FFBS algorithm, see Section 3 for numerical illustrations.

Another modification of the FFBS algorithm based on a Markov chain Monte Carlo (MCMC) sampling step was introduced in [LBS16, Section 5.2]. Instead of sampling from (12), [LBS16, Section 5.2] proposed to draw a forward path in and a sate in according to:

 ˜q(a1:i|ai+1:n,y1:n)=N∑k=1˜ϑki−1˜q(ai|ak1:i−1,ai+1:n,y1:n)δak1:i−1(a1:i−1),