Filtering with Libbi produces biased results

[Arsen Mamikonyan]

Hi all,

My goal is to use Libbi for state space modeling, and running filtering for non-linear non-gaussian models, but in order to get going, I decided to implement a simple model and see if it performs well.

So I am using a simple random-walk like Kalman filter, i.e.

model KalmanFilter {
  param q, r;
  noise eta;
  state x;
  obs y;

  sub parameter {
    q <- 0.001;
    r <- 0.001;
  }

  sub initial {
    x ~ gaussian(0, 0.001);
  }

  sub transition {
    eta ~ gaussian();
    x <- x + q * eta;
  }

  sub observation {
    y ~ gaussian(x, r);
  }
}

And I run a filtering task, on the sp500 data I found in https://github.com/lawmurray/StochasticVolatility

```

libbi filter \
  --model-file KalmanFilter.bi \
  --obs-file data/sp500.nc \
  --filter bootstrap \
  --start-time 0 \
  --end-time 100 \
  --nparticles 43210 \
  --output-file filtered.nc
```

So my concern is that the results I get are different from theoretical results [see attached] (which I calculate using a Kalman filter from filterpy python library), and libbi seems to be filtering more aggressively than the theoretical results, is this a kind of behavior I should expect?

Question 2: I tried using a Metropolis-Hastings Resampler, but it doesn't seem possible. Is it possible to use a resampler for the particles? What block should I provide to achieve this other than `--resampler` option to the script?

My code to generate the plot can be found here - https://github.com/hilearn/libbi-KalmanFilter

Cheers,

Arsen

[Dominic Steinitz]

With that number of particles, I would expect pretty close agreement between a particle filter and a Kalman filter. There are more than 3 lines and these are unlabelled. Perhaps you can remove the lines or enhance the legend?

I’d love to look at this but I am tied up today, tomorrow and some of Saturday. Perhaps you could also post your python code or maybe create a little GitHub repo so we can try to reproduce your results?

Dominic Steinitz

dom...@steinitz.org

http://idontgetoutmuch.wordpress.com

Twitter: @idontgetoutmuch

--
You received this message because you are subscribed to the Google Groups "LibBi Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to libbi-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<results.png>

[Dominic Steinitz]

Oops sorry you did post a GitHub repo

Dominic Steinitz

dom...@steinitz.org

http://idontgetoutmuch.wordpress.com

Twitter: @idontgetoutmuch

[Arsen Mamikonyan]

Unlabeled lines were +- std lines, I made them more transparent and added labels [see attached plot]

Updated the repo with the script changes. Let me know if you have trouble replicating my results.

Also, I came across this paper - http://www.matematik.lu.se/matstat/staff/jimmy/Publications/bootstrap_filter_bias.pdf

I haven't had a chance to read and understand, have you? There should not be any bias in this simple case, correct?

Cheers,

Arsen

[Dominic Steinitz]

One quick thought: are the python and libbi normal distributions parameterised the same way? Some folks use variance and others use standard deviation. That’s certainly tripped me up a few times.

Dominic Steinitz

dom...@steinitz.org

http://idontgetoutmuch.wordpress.com

Twitter: @idontgetoutmuch

<better_plot.png>

[Arsen Mamikonyan]

I've looked into that, it doesn't seem to be the issue.

I made appropriate conversions (use variance for python, std for libbi). Also the ratio of observation and process std's is what matters anyways for the Kalman filter, and that's definitely the same. I take all stds = 0.001

[Dominic Steinitz]

Hi Arsen,

I’ve managed to have a quick look after some problems building a 2 year old repo(!)

With a Kalman filter written in Haskell, I seem to get the same results as your Python code but I need to investigate further to be sure.

I have a particle filter written in Haskell so I will try that next after I have done my investigations with the Kalman filters (in Haskell and Python).

One thing that puzzles me is

  sub proposal_parameter {

    q ~ log_gaussian(q, 1);

    r ~ log_gaussian(r, 1);

  }

In the Python, these are known and not estimated. But maybe I have misunderstood.

I don’t know when I will get more time to work on this much as I would like to.

Dominic Steinitz

dom...@steinitz.org

http://idontgetoutmuch.wordpress.com

Twitter: @idontgetoutmuch

On 3 Oct 2018, at 14:55, Arsen Mamikonyan <ar...@hilearn.io> wrote:

Hi Dominic,

did you have a chance to look at the repo?

Cheers,

Arsen

[Arsen Mamikonyan]

Hi Dominic,

Thanks for looking at this.

the sub proposal_parameter block is not being used for filtering. I had it there to do parameter estimation, and forgot to remove for this demo. So ignore it, I removed from the model code.

If you can point me in the direction how to debug this inside libbi, I will try to be helpful.

Cheers,

Arsen

[Dominic Steinitz]

Hi Arsen,

I’m sorry but I don’t know how to debug this in libbi and I have zero time at the moment much as I would like to investigate this myself.

Dominic Steinitz

dom...@steinitz.org

http://idontgetoutmuch.wordpress.com

Twitter: @idontgetoutmuch

[Sebastian Funk]

1) I don't think you would should expect the two to give the same
result as you're fixing the variance of the steps and observations in
the particle filter. In your plot, notice that the biggest discrepancies
are where the Kalman filter makes a big jump. Note that LibBi also
implements a Kalman filter (with --filter kalman) which give the same
results as yours (see attached plot and R code).

2) You need to use --resampler metropolis and -C with a number greater
than 0. Does it work if you try with the latest version on github?