Bootstrap filter

Question

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$

$$x_1\sim\mathcal{N}(0,\frac{\sigma^2}{1-\phi^2})$$

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?

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$.