I am trying to implement bootstrap filter and I'm trying to understand it based on Bootstrap filter/ Particle filter algorithm(Understanding)

EDIT: Following is the example that I'm trying to solve:

Consider the stochastic volatility model

$$x_t | x_{t-1}\sim\mathcal{N}(x_t; \phi x_{t-1}, \sigma^2)$$ $$y_t | x_t\sim\mathcal{N}(y_t;0, \beta^2\exp(x_t))$$

where $x_t$ denotes the underlying latent volatility and $y_t$ the observed scaled log-returns calculated by $$y_t = 100[\log(s_t) − \log(s_{t-1})]$$

The $T = 500$ observations that we consider are log-returns during a two year period.

We would like to find an estimation for marginal filtering distribution at each time index $t = 1, ... , 𝑇$ using the bootstrap particle filter with $𝑁 = 500$ particles. I make the following assumption for the initial state $x_1$


The algorithm that I use is following:

(1) Initialization:

for i=1,..N, sample $x_1\sim\mathcal{N}(0,\frac{\sigma^2}{1-\phi^2})$

(2) Importance sampling step:

a. for i=1,.. N, sample $\tilde x_t | x_{t-1}\sim\mathcal{N}(x_t; \phi x_{t-1}, \sigma^2)$

b. evaluate $\tilde w_t^i \sim \mathcal{N}(y_t;0, \beta^2\exp(x_t))$

c. Normalize weights $w^{(i)}_{t}=\frac{\tilde{w}^{(i)}_{t}}{\sum_{i=1}^{N}\tilde{w}^{(i)}_{t}}$

(3) Resample:

a. Construct the cumulative sum of weights (CSW) by computing $c_i=c_{i-1}+w_t^i$

b. Let i = 1 and draw a starting point from the uniform distribution $u_1\sim U[0, 1/N]$

c. For j = 1, …, N, Move along the CSW by making $u_j=u_1+(j-1)/N$

d. While $u_j>c_i$ make $i=i+1$

e. Assign samples $x_t^j=x_t^i$

Is this a correct algorithm? what is marginal filtering distribution and mean of the filtering distribution?

  • $\begingroup$ Where is the question? What is the connection between a code where you generate the pairs $(x_t,y_t)$ and the particle filter? $\endgroup$
    – Xi'an
    Commented Nov 5, 2017 at 15:10
  • $\begingroup$ Dear @Xi'an , function N defines a normal distribution. I'm drawing samples from a normal distribution $x_t|x_{t-1} \tilde N(x_t;\phi x_{t-1},\sigma^2)$ and then calculate the weights based on $y_t|x_t \tilde N(y_t;0,\beta^2 \exp(x_t))$ $\endgroup$
    – Blade
    Commented Nov 5, 2017 at 15:21
  • $\begingroup$ The code is ambiguous since the same notation applies to draws (for $x_t$'s) and densities (for $y_t$'s). This seems to contradict your definition in the last part of the code. Ps_I would also advise you to write this stochastic volatility model in mathematical notations outside the Matlab code. $\endgroup$
    – Xi'an
    Commented Nov 5, 2017 at 15:26
  • $\begingroup$ @Xi'an Hope this is more clear now. $\endgroup$
    – Blade
    Commented Nov 5, 2017 at 23:51

1 Answer 1


Your question

The most common way to resample is with a multinomial distribution, resampling the same number of particles as you have. In this case the indices are sampled from a $\text{Multinomial}(N, \left[\frac{\tilde{w}_1}{\sum_i \tilde{w}_i }, \ldots, \frac{\tilde{w}_N}{\sum_i \tilde{w}_i } \right])$. In R, this is accomplished with something like s <- sample(1:N, size=N, replace=TRUE, prob=w). This probabilistically discards samples with low weights, and probabilistically duplicates samples with high weights. However there are other schemes that possess less variance, such as systematic resampling, and residual sampling.

Note that the notation in the document you cite stresses the resampling of values corresponding with the indices, instead of the indices themselves. This is fine, just note that the $i$ in $\tilde{x}^{(i)}_{0:t}$ does not correspond with the $i$ in $x_{0:t}^{(i)}$.

Some things you appear to be confused about

There are a few mistakes you seem to be making. First, I changed the notation in your question a bit. Keep in mind that you are drawing samples $X_t^i|x_{t-1} \sim N(x_t;\phi x_{t-1},\sigma^2)$ for $i=1,\ldots,N$. And for each of these samples, you revise a weight using the observation density $g(y_t|x_t^i) = N(y_t;0,\beta^2 \exp(x_t^i))$. So you're doing two things, drawing random samples, and then using those samples to evaluate a density for your weights. You use this density and evaluate it using the most recent $y_t$, and the corresponding sample $x_t^i$. At any given time you should have an $N$ dimensional array of all your weights, one for each state path, and $N$ samples that are $t$ long.

Also, I'm not sure if your code is pseudo-code, or perhaps Matlab (which I'm not familiar with), but the line y(t,i)=N(0,theta(3)^2*exp(x(t,i))); % importance weights (observation density) signifies confusion, because you should not call your weights $y_t$.

  • $\begingroup$ with current edits, let me know if you still believe that I'm confused. $\endgroup$
    – Blade
    Commented Nov 5, 2017 at 23:52
  • $\begingroup$ @Blade yeah still looks like it. I'd suggest trying to break it down into several smaller questions that are more specific. You'll get better answers, and third parties will get more out of it looking at this as a reference. $\endgroup$
    – Taylor
    Commented Nov 6, 2017 at 0:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.