ar
X
iv
:0
80
6.
42
42
v1
[
ma
th.
ST
]
26
Ju
n 2
00
8
Online data processing: comparison of Bayesian regularized
particle filters
Roberto Casarin
Department of Economics, University of Brescia
Jean-Michel Marin∗
INRIA Futurs, Projet select, Universite´ Paris-Sud
June 26, 2008
Abstract
The aim of this paper is to compare three regularized particle filters in an online data
processing context. We carry out the comparison in terms of hidden states filtering and
parameters estimation, considering a Bayesian paradigm and a univariate stochastic volatil-
ity model. We discuss the use of an improper prior distribution in the initialization of the
filtering procedure and show that the regularized Auxiliary Particle Filter (APF) outper-
forms the regularized Sequential Importance Sampling (SIS) and the regularized Sampling
Importance Resampling (SIR).
Keywords: Online data processing; Bayesian estimation; regularized particle filters; stochas-
tic volatility model
1 Introduction
The analysis of phenomena, which evolve over time is a common problem to many fields like
engineering, physics, biology, statistics, economics and finance. A time varying system can be
represented through a dynamic model, which is constituted by an observable component and an
unobservable internal state. The hidden states (or latent variables) represent the informations
we want to extrapolate from the observations.
In time series analysis, many approaches have been used for the estimation of dynamics
models. The seminal works of Kalman (1960) and Kalman and Bucy (1960) introduce filter-
ing techniques (the Kalman-Bucy filter) for continuous valued, linear and Gaussian dynamic
systems. Maybeck (1982) motivates the use of stochastic dynamic systems in engineering and
examines the estimation problems for state space models, in both a continuous and a discrete
time framework. In economics, Harvey (1989) studies the state space representation of dynamic
structural models and uses Kalman filter for hidden states filtering. Hamilton (1989) analyzes
nonlinear time series models and introduces a filter (Hamilton-Kitagawa filter) for discrete time
and discrete valued dynamic systems with a finite number of states.
In this paper, the online data processing problem is considered. In these situations, as
pointed out by Liu and Chen (1998), Markov Chain Monte Carlo (MCMC) samplers are too
much time demanding. To overcome this difficulty, some sequential Monte Carlo techniques have
∗Corresponding author: Universite´ d’Orsay, Laboratoire de Mathe´matiques (Baˆt. 425), 91405 Orsay Cedex
jean-michel.marin@inria.fr
1

been recently developed. Doucet et al. (2001) provide the state of the art on these methods.
They discuss both applications and theoretical convergence of the algorithms.
The contribution of this work is the comparison of three types of regularized particle fil-
ters - the regularized Sequential Importance Sampling (SIS), the regularized Sampling Im-
portance Resampling (SIR) and the regularized Auxiliary Particle Filter (APF) - when the
model parameters are unknown. The online estimation of model parameters is a difficult task
(Kitagawa (1998); Storvik (2002); Berzuini and Gilks (2001); Fearnhead (2002); Djuric et al.
(2002); Storvik (2002); Andrieu and Doucet (2003); Doucet and Tadic (2003); Polson et al. (2002)).
We consider here the Bayesian paradigm and the regularization (see Chen and Haykin (2002))
approach of Oudjane (2000); Liu and West (2001); Musso et al. (2001); Rossi (2004) based on a
kernel approximation in the parameter-augmented state space. We also discuss the initialization
of the filtering procedure.
This work is structured as follow. Section 2 introduces the general representation of a
Bayesian dynamic model and presents the stochastic volatility model. Section 3 reviews some
regularized particle filters. Finally, Section 4 presents the results.
2 Bayesian Dynamic Models
We introduce the general formulation of a Bayesian dynamic model and show some fundamental
relations for Bayesian inference on it. Our definition of dynamic model is general enough to
include the models analyzed in Kalman (1960), Hamilton (1994), Carter and Kohn (1994),
Harrison and West (1989) and in Doucet et al. (2001). Throughout this work, we use a notation
similar to that one commonly used in particle filters literature (see Doucet et al. (2001)).
We denote by {xt; t ∈ N}, xt ∈ X ⊆ Rnx, the hidden states of the system, by {yt; t ∈ N0},
yt ∈ Y ⊆ Rny , the observable variables and by {θt; t ∈ N}, θt ∈ Θ ⊆ Rnθ , the parameters of the
model. We denote by x0:t = (x0, . . . ,xt) the collection of hidden states up to time t and with
x−t = (x0, . . . ,xt−1,xt+1, . . . ,xT ) the collection of all hidden states without the t-th element.
We use the same notations for the observable variables and parameters.
The Bayesian state space representation of a dynamic model is given by:
yt ∼ p(yt|xt,θt,y1:t−1) measurement density ,
(xt,θt) ∼ p(xt,θt|x0:t−1,θ0:t−1,y1:t−1) transition density ,
x0 ∼ p(x0|θ0) initial density ,
θ0 ∼ pi(θ0) prior density ,
for t = 1, . . . , T .
We suppose that p(xt,θt|x0:t−1,θ0:t−1,y1:t−1) = p(xt,θt|xt−1,θt−1,y1:t−1). We also assume
that the parameters are constant over time: the transition density of the parameters is then
δθt−1(θt) with initial value θ0 = θ, δx(y) denotes the Dirac’s mass centered in x.
In that case, the joint transition of hidden states and parameters is:
p(xt,θt|xt−1,θt−1,y1:t−1) = p(xt|xt−1,θt,y1:t−1)δθt−1(θt) .
Let us denote by zt = (xt,θt) the parameter-augmented state vector and by Z the corre-
sponding augmented state space. For such models, we are interested in the prediction, filtering
2

and smoothing densities which are given by:
p(zt+1|y1:t) =
∫
Z
p(xt+1|xt,θt+1,y1:t)δθt(θt+1)p(zt|y1:t)dzt , (1)
p(yt+1|y1:t) =
∫
Z
p(yt+1|zt+1,y1:t)p(zt+1|y1:t)dzt+1 ,
p(zt+1|y1:t+1) =
p(yt+1|zt+1,y1:t)p(zt+1|y1:t)
p(yt+1|y1:t)
, (2)
p(zs|y1:t) = p(zs|y1:s)
∫
Z
p(zs+1|zs,y1:s)p(zs+1|y1:t)
p(zs+1|y1:s)
dzs+1, s < t .
Due to the high number of integrals that must be solved, previous densities may be difficult to
evaluate with general dynamics. Some Monte Carlo simulation methods, such as particle filters,
allow us to overcome these difficulties.
As an example, let us consider the stochastic volatility model. Two of the main features
of the financial time series are time varying volatility and clustering phenomena in volatility.
Stochastic volatility models widely used in finance have been introduced, in order to account
for these features. Let yt be the observable variable with time varying volatility and xt the
stochastic log-volatility process. An example of stochastic volatility model is:
yt|xt ∼ N (0, ext)
xt|xt−1,θ ∼ N
(
α+ φxt−1, σ2
)
x0|θ ∼ N
(
0, σ2/(1 − φ2)
)
θ ∼ pi(θ)
where θ = (α, log((1 + φ)/(1− φ)), log(σ2)). The choice of pi(θ) will be discussed in Section 4.
Fig. 1 shows two simulated paths of yt and xt.
In the next section, we deal with the problem of parameters and states joint estimation in
a kernel-regularized sequential Monte Carlo framework.
3 Regularized particular filters
For making inference on the Bayesian dynamic model given in Section 2 in an online data
processing context, MCMC algorithms are too much time demanding. Sequential impor-
tance sampling and more advanced sequential Monte Carlo algorithms called Particle Filters
(Doucet et al., 2001) represent a promising alternative. The main advantage in using particle
filters is that they can deal with nonlinear models and non-Gaussian innovations. In contrast
to Hidden Markov Model filters, which work on a state space discretized to a fixed grid, par-
ticle filters focus sequentially on the higher density regions of the state space. This feature is
common to one of the early sequential methods, the Adaptive Importance Sampling algorithm
due to West (1992, 1993).
Different particle filters exist in the literature and different simulation approaches like rejec-
tion sampling, MCMC and importance sampling, can be used for the construction of a particle
filter. In this work, we present some kernel-regularized particle filters, which combine the im-
portance sampling reasoning with a suitable modification of the importance weights. The reg-
ularization approach we use is the same than the one of Liu and West (2001) and Musso et al.
(2001). This approach relies upon a kernel-based reconstruction of the empirical filtering den-
sities which produces a systematic modification of the true importance weights.
3

0 200 400 600 800 1000
−
6
−
2
0
2
4
Daily
Time
0 200 400 600 800 1000
−
6
−
2
0
2
4
Weekly
Time
Figure 1: Simulation paths for xt (grey line) and yt (black line). Upper plot: daily dataset
(α = 0, φ = 0.99 and σ2 = 0.01). Bottom plot: weekly dataset (α = 0, φ = 0.9 and σ2 = 0.1).
3.1 Regularized SIS
Let us start from the non-regularized SIS. We assume that at iteration t > 0 a properly weighted
particle set {xit,θit, γit}Ni=1, approximating the filtering density p(xt,θt|y1:t), is available. The
empirical distribution corresponding to this approximation is:
pN (xt,θt|y1:t) =
N
∑
i=1
γitδ(xit,θit)(xt,θt) . (3)
The particles set, {xit,θit, γit}Ni=1, can be viewed as a random discretisation of the state space
X ×Θ with associated probability weights {γit}Ni=1. Thanks to this discretisation, it is possible
to approximate the prediction and filtering densities given in (1) and (2):
pN (xt+1,θt+1|y1:t) =
N
∑
i=1
γitp(xt+1|xit,θt+1,y1:t)δθit(θt+1) ,
pN(xt+1,θt+1|y1:t+1) ∝
N
∑
i=1
γitp(yt+1|xt+1,θt+1,y1:t)p(xt+1|xit,θt+1,y1:t)δθit(θt+1) .
The goal is now to obtain N particles {xit+1,θit+1, γit+1}Ni=1 from the filtering density in (2).
It is proposed to sample (xit+1,θit+1) according to the importance density q(·|xit,θit,y1:t+1). The
importance weight of particle (xit+1,θit+1) is then calculated using the recursive formula:
γit+1 ∝ γit
p(yt+1|xit+1,θit+1,y1:t)p(xit+1|xit,θit+1,y1:t)δθit(θ
i
t+1)
q(xit+1,θit+1|xit,θit,y1:t+1)
. (4)
4

The choice of an optimal importance density q(·|xit,θit,y1:t+1), that is, a density which min-
imizes the variance of the importance weights is discussed in Pitt and Shephard (1999) and
Crisan and Doucet (2000). In many cases, it is not possible to use this optimal importance
density as the weight updating associated to the this density does not admit a closed-form
expression. In that case, the transition density of the parameter-augmented state vector repre-
sents a natural alternative for the importance density. Indeed, the transition density represents
a sort of prior at time t for the parameter-augmented state vector (xit+1,θit+1).
In our case, due to the presence of the Dirac point mass at the numerator of the weights it
is impossible to modify over the filtering iterations the particle values for the parameters. In
practice due to the loss of particle diversity in the parameter space, the weights will tend to
zeros and of course stay zero for ever, so we are facing a problem of degeneracy of the empirical
filtering distribution. This scenario motivates particle filtering methods known as regularized
particle filters. In order to avoid the degeneracy problem and to force the exploration of the
parameter space toward regions which are not covered by the prior distribution, Liu and West
(2001) and Musso et al. (2001) propose to use a regularized version of the filtering density. This
approach results in the modification of the weights in (4) and the definition of a new set of
weights:
ωit+1 ∝ ωit
p(yt+1|xit+1,θit+1,y1:t)p(xit+1|xit,θit+1,y1:t)Kh(θit+1 − θit)
q(xit+1,θit+1|xit,θit,y1:t+1)
where Kh(y) = h−dK(y/h) is a regularization kernel, K being a positive function defined on
Rnθ and h a positive smoothing factor (bandwidth).
The modification of the importance weights defined in (4) results from two steps. The first
one is the regularization of the empirical density in (3) by a kernel estimator:
pRN (xt,θt|y1:t) =
N
∑
i=1
ωitδxit(xt)Kh(θt − θ
i
t) .
The second one is the application of an importance sampling argument to the approximated
filtering density:
pRN (xt+1,θt+1|y1:t+1) =
N
∑
i=1
ωitp(yt+1|xt+1,θt+1,y1:t)p(xt+1|xit,θt+1,y1:t)Kh(θt+1 − θit) .
The convergences results associated with this type of approximation are recalled in Musso et al.
(2001) and Oudjane (2000). Under some usual conditions on the kernel, when the number of
particles increases to infinity, the regularized empirical density converges to the right one for
various criteria. For instance, we have pRN −→L2 p.
Thanks to this approximation, the regularization kernel becomes the natural choice for the
parameters proposal distribution. Thus, we sample (xit+1,θit+1) according to:
q(xt+1|xit,θt+1,y1:t+1)Kh(θt+1 − θit) .
In that case, we have:
ωit+1 ∝ ωit
p(yt+1|xit+1,θit+1,y1:t)p(xit+1|xit,θit+1,y1:t)
q(xit+1|xit,θit+1,y1:t+1)
.
In Algorithm 1, we give a pseudo-code representation of this method.
5

Algorithm 1 - Regularized SIS Particle Filter -
· At time t = 0, for i = 1, . . . , N , simulate zi0 ∼ p(z0) and set ωi0 = 1/N
· At time t > 0, given {xit,θit, ωit}Ni=1, for i = 1, . . . , N :
1. Simulate θit+1 ∼ Kh(θt+1 − θit)
2. Simulate xit+1 ∼ q(xt+1|xit,θit+1,y1:t+1)
3. Update the weights: ωit+1 ∝ ωit
p(yt+1|xit+1,θit+1,y1:t)p(xit+1|xit,θit+1,y1:t)
q(xit+1|xit,θit+1,y1:t+1)
.
3.2 Regularized SIR
As it is well known in the literature (see for example Arulampalam et al. (2001)), basic SIS
algorithms have a degeneracy problem. After some iterations the empirical distribution degen-
erates into a Dirac’s mass on a single particle. This due to the fact that the variance of the
importance weights is non-decreasing over time (see Doucet et al. (2000)). In order to solve this
degeneracy problem, Gordon et al. (1993) introduce the SIR algorithm. This algorithm belongs
to a wider class of bootstrap filters. At each iteration, a resampling step is used to generate
a new set of particles. After this resampling step, the weights of the resampled particles are
uniformly distributed over the particle set.
In the initial SIR, the resampling step is done at each iteration of the algorithm. This
systematic resampling can introduce extra Monte Carlo variations, see Liu and Chen (1998).
This can be reduced be doing resampling only when the Effective Sample Size (ESS) is small.
The ESS measures the overall efficiency of an importance sampling algorithm. The ESS is a
function of the coefficient of variation of the importance weights. At iteration t, the empirical
EES is
ESSt =
N
1 +N
N
∑
i=1
(
ωit −N−1
N
∑
i=1
ωit
)2
/
(
N
∑
i=1
ωit
)2 .
In Algorithm 2, we give a pseudo-code representation of this method.
Algorithm 2 - Regularized SIR Particle Filter -
· At time t = 0, for i = 1, . . . , N , simulate zi0 ∼ p(z0) and set ωi0 = 1/N
· At time t > 0, given {xit,θit, ωit}Ni=1, for i = 1, . . . , N :
1. Simulate θit+1 ∼ Kh(θt+1 − θit)
2. Simulate xit+1 ∼ q(xt+1|xit, θ˜
i
t+1,y1:t+1)
3. Update the weights: ωit+1 ∝ ωti
p(yt+1|x˜it+1, θ˜
i
t+1,y1:t)p(x˜
i
t+1|xit, θ˜
i
t+1,y1:t)
q(x˜it+1|xit, θ˜
i
t+1,y1:t+1)
4. If ESSt+1 < κ, simulate {xit+1,θit+1}Ni=1 from {xit+1,θit+1, ωit+1}Ni=1 (Multi-
nomial resampling) and set ωit+1 = 1/N .
The value of κ < N is calibrated depending on the problem.
6

3.3 Regularized APF
Due to the resampling step, the basic SIR algorithm produces a progressive impoverishment (loss
of diversity) of the information contained in the particle set. To overcome this difficulty, many
solutions have been proposed in the literature. We refer to the APF due to Pitt and Shephard
(1999) and to the regularized APF algorithm due to Liu and West (2001). In order to avoid
the resampling step, the APFs use the particle index (auxiliary variable) to select most rep-
resentative particles in the proposal of the new particles. The regularized joint distribution of
parameter-augmented state vector and the particle index is:
pRN (xt+1,θt+1, i|y1:t+1) ∝ p(yt+1|xt+1,θt+1,y1:t)p(xt+1|xit,θit,y1:t)Kh(θt+1 − θit)ωit .
A sample approximating that distribution can be obtained by using the proposal:
q(xit+1,θit+1, ji|y1:t+1) = p(xit+1|xj
i
t ,θit+1,y1:t)Kh(θit+1 − θj
i
t )q(j
i|y1:t+1)
where
q(ji|y1:t+1) ∝ p(yt+1|µj
i
t+1,m
ji
t+1,y1:t)w
ji
t ,
µj
i
t+1 and m
ji
t+1 are evaluated using the initial particle set. Therefore, the importance weight of
particle (xit+1,θit+1, ji) is:
ωit+1 ∝
p(yt+1|xit+1,θit+1,y1:t)
p(yt+1|µj
i
t+1,m
ji
t+1,y1:t)
.
In Algorithm 3 we give a pseudo-code representation of the regularized APF.
Algorithm 3 - Regularized Auxiliary Particle Filter -
· At time t = 0, for i = 1, . . . , N , simulate zi0 ∼ p(z0) and set ωi0 = 1/N
· At time t > 0, given {xit,θit, ωit}Ni=1, for i = 1, . . . , N :
1. Simulate ji ∼ q(j|y1:t+1) with j ∈ {1, . . . , N} (Multinomial sampling)
where µjt+1 = E(xt+1|xjt ,θjt ) and mjt+1 = E(θt+1|θjt )
2. Simulate θit+1 ∼ Kh(θt+1 − θj
i
t )
3. Simulate xit+1 ∼ p(xt+1|xj
i
t ,θit+1,y1:t)
4. Update particles weights: ωit+1 ∝
p(yt+1|xit+1,θit+1,y1:t)
p(yt+1|µj
i
t+1,m
ji
t+1,y1:t)
.
We can say that, in the APF, the selection step is done before simulating the hidden states.
This selection depends on the current value of the observable. Therefore,
• the APF is a standard way to construct a proposal distribution for the hidden states that
depends on the current of the particle;
• as we will see after, to use this selection step and the transition distribution as proposal
distribution for the hidden states, results in a good proposal distribution.
In the next section, we compare the performances of some regularized SIS, SIR and APF
for the stochastic volatility model.
7

4 Application to the stochastic volatility model and comparison
In this section, we apply the three regularized particle filters to the stochastic volatility model
presented in Section 2. We assume that the initial value of the SV process follows the stationary
distribution:
x0 ∼ N (0, σ2/(1 − φ2)) .
For the parameters β, σ and φ, we assume the prior
p(β2, φ, σ2) = 1/(σβ)I(−1,1)(φ) ,
where β = eα. We constrain the parameter φ to take values in the open interval (−1, 1) in
order to impose the usual stationarity condition. As we use an improper prior, it is not possible
to use the prior distribution for initializing the three particle filters. We need to start with a
proper weighted sample {xi0,θi0, ωi0}Ni=1. We propose to:
1) start the sequential filtering procedure at least at the value n of t, such that the posterior
distribution of the parameters given all the observations up to time n is well defined: for
the considered SV model, this corresponds to set n ≥ 2;
2) use a Markov Chain Monte Carlo (MCMC) algorithm to create a sample with uniform
weights.
For the prior given above, the full conditional distributions are
β2| · · · ∼ IG
(
n
∑
t=1
y2t exp(−xt)/2, (n − 1)/2
)
,
σ2| · · · ∼ IG
(
n
∑
t=2
(xt − φxt−1)2/2 + x21(1− φ2), (n − 1)/2
)
,
pi(φ| · · · ) ∝ (1− φ2)1/2 exp
(
−φ2
n−1
∑
t=2
x2t − 2φ
n
∑
t=2
xtxt−1
)
/2σ2I(−1,+1)(φ) ,
pi(xt| · · · ) ∝ exp
{
− 1
2σ2
(
(xt − α− φxt−1)2 − (xt+1 − α− φxt)2
)
− 1
2
(
xt + y2t exp(−xt)
)
}
.
The full conditional distributions of φ and xt are not conventional and the standard Gibbs
does not apply. We propose to use the Metropolis-Hastings within Gibbs algorithm studied
in Celeux et al. (2006). A detailed description of the proposal distributions for φ and xt can
be found in Celeux et al. (2006). In that paper, the authors compare this MCMC scheme
to an iterated importance sampling one. Note that one could alternatively use this iterated
importance sampling algorithm to create a first weighted sample.
Given the initial weighted random sample
{
xit,θit, ωit
}N
i=1, where θt = (αt, log((1 + φt)/(1−
φt)), log(σ2t )), if we use the transition density as proposal distribution for the hidden states, the
regularized SIS performs the following steps:
8

200 400 600 800
−
1.
5
0.
0
1.
5 SIS
200 400 600 800
−
2
0
2
SIS
200 400 600 800
−
1.
5
0.
0
1.
5 SIR−r−p
200 400 600 800
−
2
0
2
SIR−r−p
200 400 600 800
−
1.
5
0.
0
1.
5 APF
200 400 600 800
−
2
0
2
APF
Figure 2: Daily (left column) and weekly (right column) true (black line) and filtered (grey line)
log-volatility.
For n ≤ t ≤ T − 1 and for i = 1, . . . , N:
(i) Simulate θit+1 ∼ N
(
aθit + (1− a)θ¯t, b2Vt
)
where Vt and θ¯t are the
empirical covariance matrix and the empirical mean respectively,
a ∈ [0, 1] and b2 = (1− a2),
(ii) Simulate xit+1 ∼ N
(
αit+1 + φ
i
t+1x
i
t,
(
σ2
)i
t+1
)
,
(iii) Update the weights as follow
wit+1 ∝ wit exp
{
−1
2
[
y2t+1 exp
(
−xit+1
)
+ xit+1
]
}
.
In the following, we call the previous scheme SIS. Instead of the transition density, we can use
a proposal distribution which depends on the current value of the observable. For instance,
we can resort to the proposal distribution used for the hidden states, in the MCMC initializa-
tion step. This proposal has been introduced by Shephard and Pitt (1997). It is based on a
quadratic Taylor expansion of exp(xt), more details can be found in Shephard and Pitt (1997)
or Celeux et al. (2006). It defines a new SIS which is called SIS-p in the following.
The results of a typical run of the regularized SIS on the synthetic dataset in Fig. 1, with
N = 10, 000 particles and n = 100 for the Gibbs initialization, are given in Fig. from 2 to 6. We
can see (last row in Fig. 2) that after a few iterations the filtered log-volatility does not fit well
to the true log-volatility. We measure sequentially the filtering performance of the regularized
SIS by evaluating the cumulated Root Mean Square Error (RMSE). It measures the distance
9

between the true and the filtered states and is defined as: RMSEt = {1t
∑t
u=1(zˆu − zu)2}
1
2 ,
where zˆt is the filtered state, which includes also the parameter sequential estimate. The RMSEs
cumulate rapidly over time in both daily and weekly datasets (see upper and bottom plots in
Fig. 3). The poor performance of the regularized SIS is due to the fact that the empirical
posterior of the states and parameters degenerates into a Dirac’s mass after a few iterations.
The ESSs in Fig. 4 show that the regularized SIS degenerates after 30 iterations in both the
daily and weekly datasets. We give some results on SIS-p in the following.
If we use the transition density as proposal distribution for the hidden states, the regularized
SIR performs the following step:
For n ≤ t ≤ T − 1 and for i = 1, . . . , N:
(i) Simulate θit+1∼N
(
aθit + (1− a)θ¯t, b2Vt
)
where Vt and θ¯t are the
empirical covariance matrix and the empirical mean respectively
and a ∈ [0, 1] and b2 = (1− a2),
(ii) Simulate xit+1 ∼ N
(
αit+1 + φ
i
t+1x
i
t,
(
σ2
)i
t+1
)
,
(iii) Update the weights
wit+1 ∝ wit exp
{
−1
2
[
y2t+1 exp{−xit+1}+ xit+1
]
}
,
(v) If ESSt+1 < κ, simulate zit+1 ∼
∑N
j=1 w
j
t+1δzjt+1
(zt+1) and set wit+1 =
1/N.
If κ = N , the resampling step is done all the time. In that case, we call SIR the previous scheme.
After some numerical experiments, we have found that a good value for κ is κ = 0.9 × N . In
that case, the resampling step is done at regular time intervals and we called SIR-r the resulting
algorithm. Moreover, as for the SIS, we can resort to the proposal of Shephard and Pitt (1997)
for the hidden states. In that case, with κ = N , the corresponding algorithm called SIR-p.
With κ = 0.9 ×N the corresponding algorithm is called SIR-r-p.
The regularized APF performs the following steps:
For n ≤ t ≤ T − 1 and for i = 1, . . . , N:
(i) Simulate ji ∼ q(j) ∝∑Nk=1 wktN (yt+1|µkt+1)δk(j) where µkt+1 = φkt xkt + αkt ,
(ii) Simulate θit+1∼N
(
aθj
i
t + (1− a)θ¯t, b2Vt
)
where Vt and θ¯t are the
empirical variance matrix and the empirical mean respectively
and a ∈ [0, 1] and b2 = (1− a2),
(iii) Simulate xit+1 ∼ N
(
xt+1|αit+1 + φit+1xj
i
t ,
(
σ2
)i
t+1
)
,
(iv) Update the weights
wit+1 ∝ exp
{
−1
2
[
y2t+1
(
exp{−xit+1} − exp{−µj
i
t+1}
)
+ xit+1 − µj
i
t+1
]
}
.
10

200 400 600 800 1000
0.
0
0.
4
0.
8
1.
2
Daily
200 400 600 800 1000
0.
0
0.
4
0.
8
1.
2
Weekly
Figure 3: Daily (upper plot) and weekly (bottom plot) Root Mean Square Errors for the regu-
larized APF (solid line), SIR-r-p (dashed line) and SIS (dotted line) over iterations.
200 400 600 800
0
40
00
10
00
0 SIS
200 400 600 800
0
40
00
10
00
0 SIS
200 400 600 800
0
40
00
10
00
0 SIR−r−p
200 400 600 800
0
40
00
10
00
0 SIR−r−p
200 400 600 800
0
40
00
10
00
0 APF
200 400 600 800
0
40
00
10
00
0 APF
Figure 4: Daily (left column) and weekly (right column) Effective Sample Sizes over iterations.
11

APF for α APF for φ APF for σ2
SIR-r-p for α SIR-r-p for φ SIR-r-p for σ2
Figure 5: Evolution on daily dataset of the empirical posterior distributions of α, φ and σ2.
APF for α APF for φ APF for σ2
SIR-r-p for α SIR-r-p for φ SIR-r-p σ2
Figure 6: Evolution on weekly dataset of the empirical posterior distributions of α, φ and σ2.
12

Daily Data
θ SIS SIS-p SIR SIR-p SIR-r SIR-r-p APF
α 0.00719 0.00945 0.00885 0.00925 0.00315 0.00912 0.00065
φ 0.66767 0.83264 0.12433 0.13355 0.15252 0.13456 0.00855
σ2 0.89327 0.87910 0.00676 0.00670 0.00643 0.00654 0.00506
Table 1: Mean Square Errors of the estimators of α, φ and σ2. The Mean Square Errors are
estimated using the last iteration of the 10 independent runs of the filters.
Weekly Data
θ SIS SIS-p SIR SIR-p SIR-r SIR-r-p APF
α 0.00534 0.00487 0.00589 0.00442 0.00380 0.00431 0.00016
φ 0.51290 0.55648 0.05292 0.03754 0.04885 0.03687 0.00029
σ2 0.70540 0.71242 0.00010 0.00009 0.00009 0.00009 0.00008
Table 2: Mean Square Errors of the estimators of α, φ and σ2. The Mean Square Errors are
estimated using the last iteration of the 10 independent runs of the filters.
Note that, following Pitt and Shephard (1999), one could alternatively use in the selection step
a value of µkt+1 based on the Taylor expansion of the likelihood at time t + 1. For the three
regularized particle filters, we have used a Gaussian kernel where the parameter a is fixed
following the usual optimal criterion.
We apply the regularized SIR-r-p and APF with N = 10, 000 and n = 100 to the weekly
and daily datasets of Fig. 1 and obtain the results given in Fig. from 2 to 6. The regularized
SIR-r-p and APF outperform the regularized SIS in terms of ESSs and cumulated RMSEs.
The ESSs can detect the degeneracy in the particle weights, but is not useful to determine the
presence of another form of degeneracy, that is the absence of diversity in the particle values.
The histogram of the empirical filtering distribution allows us to detect this second form of
degeneracy.
As our work deals with the sequential estimation of the parameters, we choose to show
the histogram of the parameters posterior. Fig. 5 and 6 exhibit the evolution over the filters
iterations of the posterior of the parameters α, φ and σ2. In both the daily and the weekly
datasets, after a few iterations the empirical posterior of the regularized SIR-r-p degenerates
into a Dirac’s mass.
To confirm the previous results, we have done ten independent runs of the seven algorithms:
SIS, SIS-p, SIR, SIR-p, SIR-r, SIR-r-p and APF. The weekly and daily datasets vary across
the 10 experiments. Our simulation study confirms the results of the single-run experiment.
Fig. 7 and 8 show a comparison between the regularized schemes in terms of RMSEs. The
RMSEs are estimated over the 10 independent runs of the algorithms. The regularized SIR-r-p
and APF outperform the others algorithms in both the daily and the weekly datasets. The
estimated Mean Square Errors for the parameters α, φ and σ2 (see Table 1 and 2), based on
10 independent runs of the filters, show that the regularized APF outperforms all the others
schemes in term of parameters estimation.
Fig. 9 and 10 show a comparison of the filters in terms of ESS. As one could expect, in all the
independent runs the regularized SIS and SIS-p weights degenerate after a few iterations. We
also observe that the ESSs of the regularized SIR, SIR-r and APF are substantially equivalent.
On the other hand, the average ESSs of the SIR-p and SIR-r-p are lower than the ones of the
13

SIR, SIR-r and APF.
5 Conclusion
In this work we bring into action the kernel regularization technique for particle filters and deal
with the online parameter estimation problem. While the regularized APF has been already
used for the parameter estimation, the regularized versions of SIS and SIR have not been
considered to that aim. We focus on the joint estimation of the states and parameters and
compare some algorithms on a Bayesian nonlinear model: the Bayesian stochastic volatility
model. As we expected, we find evidence of the degeneracy of two different regularized SIS.
Finally, we find that, in terms of parameters estimation, the regularized APF outperforms all
the others schemes.
Acknowledgments
We thank Professor Christian Robert for its helpful comments and a careful reading of a pre-
liminary version of the paper.
References
Andrieu, C. and Doucet, A. (2003). Online expectation–maximization type algorithms for
parameter estimation in general state space models. Proc. IEEE ICASSP, 6:VI–69–VI–72.
Arulampalam, S., Maskell, S., Gordon, N., and Clapp, T. (2001). A tutorial on particle filters
for on-line nonlinear/non-gaussian bayesian tracking. Technical Report, QinetiQ Ltd., DSTO,
Cambridge.
Berzuini, C. and Gilks, W. (2001). Following a moving average Monte Carlo inference for
dynamic Bayesian models. Journal of Royal Statistical Society, B, 63:127–146.
Carter, C. and Kohn, R. (1994). On Gibbs Sampling for State Space Models. Biometrika,
81(3):541–553.
Celeux, G., Marin, J.-M., and Robert, C. (2006). Iterated importance sampling in missing data
problems. Computational Statistics and Data Analysis, 50(12):3386–3404.
Chen, Z. and Haykin, S. (2002). On different facets of regularization theory. Neural Comput.,
14:2791–2846.
Crisan, D. and Doucet, A. (2000). Convergence of sequential monte carlo methods. Technical
Report N. 381, CUED-F-INFENG.
Djuric, P. M., Kotecha, J., Esteve, F., and Perret, E. (2002). Sequential parameter estimation
of time-varying non-gaussian autoregressive processes. EURASIP Journal on Applied Signal
Processing, 8:865–875.
Doucet, A., de Freitas, N., and Gordon, N. (2001). Sequential Monte Carlo Methods in Practice.
Springer-Verlag.
Doucet, A., Godsill, S., and Andrieu, C. (2000). On sequential Monte Carlo sampling methods
for Bayesian filtering. Statistics and Computing, 10:197–208.
14

Doucet, A. and Tadic, C. (2003). Parameter estimation in general state-space models using
particle methods. Annals of the Institute of Statistical Mathematics, 55(2):409–422.
Fearnhead, P. (2002). MCMC, sufficient statistics and particle filter. Journal of Computational
and Graphical Statistics, 11:848–862.
Gordon, N., Salmond, D., and Smith, A. F. M. (1993). Novel approach to nonlinear and non-
gaussian bayesian state estimation. IEE Proceedings-F, 140:107–113.
Hamilton, J. (1994). Time Series Analysis. Princeton University Press.
Hamilton, J. D. (1989). A new approach to the economic analysis of nonstationary time series
and the business cycle. Econometrica, 57:357–384.
Harrison, J. and West, M. (1989). Bayesian Forecasting and Dynamic Models. Springer-Verlag,
2 edition.
Harvey, A. (1989). Forecasting, structural time series models and the Kalman filter. Cambridge
University Press.
Kalman, R. (1960). A new approach to linear filtering and prediction problems. Transaction of
the ASME, Journal of Basic Engineering, Series D, 82:35–45.
Kalman, R. and Bucy, R. (1960). New results in linear filtering and prediction problems.
Transaction of the ASME, Journal of Basic Engineering, Series D, 83:95–108.
Kitagawa, G. (1998). A self-organized state-space model. Journal of the American Statistical
Association, 93(443):1203–1215.
Liu, J. and Chen, R. (1998). Sequential Monte Carlo methods for dynamical system. Journal
of the American Statistical Association, 93:1032–1044.
Liu, J. and West, M. (2001). Combined parameter and state estimation in simulation based
filtering. In Doucet, A., de Freitas, N., and Gordon, N., editors, Sequential Monte Carlo
Methods in Practice. Springer-Verlag.
Maybeck, P. (1982). Stochastic Models, Estimation and Control, volume 1-3. Academic Press.
Musso, C., Oudjane, N., and Legland, F. (2001). Improving regularised particle filters. In
Doucet, A., de Freitas, N., and Gordon, N., editors, Sequential Monte Carlo Methods in
Practice. Springer-Verlag.
Oudjane, N. (2000). Stabilite´ et approximation particulaires en filtrage non-line´aire. Application
au pistage. The`se du Doctorat en Science, Universite´ de Rennes.
Pitt, M. and Shephard, N. (1999). Filtering via Simulation: Auxiliary Particle Filters. Journal
of the American Statistical Association, 94(446):590–599.
Polson, N. G., Stroud, J. R., and Mu¨ller, P. (2002). Practical Filtering with sequential parameter
learning. Tech. report, Graduate School of Business, University of Chicago.
Rossi, V. (2004). Filtrage non line´aire par noyaux de convolution. Application a` un proce´de´
de de´pollution biologique. The`se du Doctorat en Science, E´cole Nationale Supe´rieure
Agronomique de Montpellier.
15

Shephard, N. and Pitt, M. (1997). Likelihood Analysis of Non-Gaussian Measurement Time
Series. Biometrika, 84:653–667.
Storvik, G. (2002). Particle filters for state space models with the presence of unknown static
parameters. IEEE Transactions on Signal Processing, 50:281–289.
West, M. (1992). Mixture models, Monte Carlo, Bayesian updating and dynamic models.
Computer Science and Statistics, 24:325–333.
West, M. (1993). Approximating posterior distribution by mixtures. Journal of Royal Statistical
Society, B, 55:409–442.
16

200 400 600 800 1000
0.
0
1.
5 SIS
200 400 600 800 1000
0.
0
1.
5 SIS−p
200 400 600 800 1000
0.
0
1.
5 SIR
200 400 600 800 1000
0.
0
1.
5 SIR−p
200 400 600 800 1000
0.
0
1.
5 SIR−r
200 400 600 800 1000
0.
0
1.
5 SIR−r−p
200 400 600 800 1000
0.
0
1.
5 APF
Figure 7: Comparison on daily datasets of average cumulative Root Mean Square Errors between
the true and the filtered log-volatility (black line). We represent the area between maximum
and minimum cumulative Root Mean Square Errors (grey area).
17

200 400 600 800 1000
0.
0
1.
2
SIS
200 400 600 800 1000
0.
0
1.
2
SIS−p
200 400 600 800 1000
0.
0
1.
2
SIR
200 400 600 800 1000
0.
0
1.
2
SIR−p
200 400 600 800 1000
0.
0
1.
2
SIR−r
200 400 600 800 1000
0.
0
1.
2
SIR−r−p
200 400 600 800 1000
0.
0
1.
2
APF
Figure 8: Comparison on weekly datasets of average cumulative RMSEs between the true and
the filtered log-volatility (black line). We represent the area between maximum and minimum
cumulative Root Mean Square Errors (grey area).
18

200 400 600 800 1000
0
80
00
SIS
200 400 600 800 1000
0
80
00
SIS−p
200 400 600 800 1000
0
80
00
SIR
200 400 600 800 1000
0
80
00
SIR−p
200 400 600 800 1000
0
80
00
SIR−r
200 400 600 800 1000
0
80
00
SIR−r−p
200 400 600 800 1000
0
80
00
APF
Figure 9: Comparison on daily datasets of average Effective Sample Sizes (black line). We
represent the area between maximum and minimum Effective Sample Sizes (grey area).
19

200 400 600 800 1000
0
80
00
SIS
200 400 600 800 1000
0
80
00
SIS−p
200 400 600 800 1000
0
80
00
SIR
200 400 600 800 1000
0
80
00
SIR−p
200 400 600 800 1000
0
80
00
SIR−r
200 400 600 800 1000
0
80
00
SIR−r−p
200 400 600 800 1000
0
80
00
APF
Figure 10: Comparison on weekly datasets of average Effective Sample Sizes (black line). We
represent the area between maximum and minimum Effective Sample Size (grey area).
20