[pystatsmodels] State space simulation smoothing

[Alastair Heggie]

I am having trouble getting simulation smoothing of state space models to
work.

Suppose I have an UnobservedComponents model.

I can then do:

fitted = mod.fit()
sim = mod.simulation_smoother()

sim.simulate(initial_state_variates=fitted.smoothed_state[:,0],
pretransformed_variates=True)

However the resulting output does not seem to be conditional on the
observations. Is there a kwarg I need to pass to the simulation_smoother to
get this to work.

I have tried to figure this out from studying the source but to no avail.

Best wishes,
Alastair

[Chad Fulton]

Post by Alastair Heggie
I am having trouble getting simulation smoothing of state space models to
work.
Suppose I have an UnobservedComponents model.
fitted = mod.fit()
sim = mod.simulation_smoother()
sim.simulate(initial_state_variates=fitted.smoothed_state[:,0],
pretransformed_variates=True)
However the resulting output does not seem to be conditional on the
observations. Is there a kwarg I need to pass to the simulation_smoother to
get this to work.
I have tried to figure this out from studying the source but to no avail.
Best wishes,
Alastair

There are no kwargs that you should have to pass.

It's hard to say more based on your e-mail, though. Can you describe what
you're trying to do in more detail, e.g. what you're expecting and what's
actually happening, and the version of Statsmodels that you're using?

[Alastair Heggie]

Hi Chad, I'll try to be a bit more clear:

I'm trying to draw "samples of state or disturbance vectors conditional on
the observations held fixed" (quoting Durbin and Koopman).

Ultimately the use case I have in mind is a Bayesian analysis. Using pymc3
I will take draws from the posterior distribution of the error covariance
parameters and then, conditional on these parameters and the observations,
I wish to simulate the evolution of the latent state parameters.

Ignoring the Bayesian application here's some reproducible code:
%matplotlib inline
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

# Get some data
nile = sm.datasets.nile.load_pandas().data['volume']
nile.index = pd.DatetimeIndex(start="1871",end="1971", freq='y')
# Estimate the parameters of the model
mod = sm.tsa.UnobservedComponents(nile, 'llevel')
res = mod.fit()
# Get the simulation smoother and simulate
sim = res.model.simulation_smoother()
sim.simulate()
# Plot the filtered state and the state produced by the simulation smoother
pd.Series(sim.generated_obs[0,:]).plot(observed=False)
res.plot_components()

I would have expected state produced by the simulation smoother to look
similar to the estimated state from the fitted model. However, as far as I
can see this is just simulating a local level process with the given
parameters.

I have built statsmodels from the version on github (0.10.0 according to
the output of "conda list")

Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I am having trouble getting simulation smoothing of state space models to
work.
Suppose I have an UnobservedComponents model.
fitted = mod.fit()
sim = mod.simulation_smoother()
sim.simulate(initial_state_variates=fitted.smoothed_state[:,0],
pretransformed_variates=True)
However the resulting output does not seem to be conditional on the
observations. Is there a kwarg I need to pass to the simulation_smoother to
get this to work.
I have tried to figure this out from studying the source but to no avail.
Best wishes,
Alastair

There are no kwargs that you should have to pass.
It's hard to say more based on your e-mail, though. Can you describe what
you're trying to do in more detail, e.g. what you're expecting and what's
actually happening, and the version of Statsmodels that you're using?

[Chad Fulton]

Post by Alastair Heggie
I'm trying to draw "samples of state or disturbance vectors conditional on
the observations held fixed" (quoting Durbin and Koopman).
Ultimately the use case I have in mind is a Bayesian analysis. Using pymc3
I will take draws from the posterior distribution of the error covariance
parameters and then, conditional on these parameters and the observations,
I wish to simulate the evolution of the latent state parameters.
%matplotlib inline
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
# Get some data
nile = sm.datasets.nile.load_pandas().data['volume']
nile.index = pd.DatetimeIndex(start="1871",end="1971", freq='y')
# Estimate the parameters of the model
mod = sm.tsa.UnobservedComponents(nile, 'llevel')
res = mod.fit()
# Get the simulation smoother and simulate
sim = res.model.simulation_smoother()
sim.simulate()
# Plot the filtered state and the state produced by the simulation smoother
pd.Series(sim.generated_obs[0,:]).plot(observed=False)
res.plot_components()
I would have expected state produced by the simulation smoother to look
similar to the estimated state from the fitted model. However, as far as I
can see this is just simulating a local level process with the given
parameters.
I have built statsmodels from the version on github (0.10.0 according to
the output of "conda list")
Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I am having trouble getting simulation smoothing of state space models
to work.
Suppose I have an UnobservedComponents model.
fitted = mod.fit()
sim = mod.simulation_smoother()
sim.simulate(initial_state_variates=fitted.smoothed_state[:,0],
pretransformed_variates=True)
However the resulting output does not seem to be conditional on the
observations. Is there a kwarg I need to pass to the simulation_smoother to
get this to work.
I have tried to figure this out from studying the source but to no avail.
Best wishes,
Alastair

There are no kwargs that you should have to pass.
It's hard to say more based on your e-mail, though. Can you describe what
you're trying to do in more detail, e.g. what you're expecting and what's
actually happening, and the version of Statsmodels that you're using?

Hi Alastair,

Thanks for posting that!

I think maybe you're not looking at the variable you want - I'm guessing
that you're looking for `sim.simulated_state`, which is the draw from the
joint posterior of the states.

Best,
Chad

[Alastair Heggie]

Great, that works, not sure how I missed that...

For reference here's a notebook summarizing the use case I had in mind

https://gist.github.com/AliHeggie/3f2f56963a0e6893043d4116dfa3caed

Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I'm trying to draw "samples of state or disturbance vectors conditional
on the observations held fixed" (quoting Durbin and Koopman).
Ultimately the use case I have in mind is a Bayesian analysis. Using
pymc3 I will take draws from the posterior distribution of the error
covariance parameters and then, conditional on these parameters and the
observations, I wish to simulate the evolution of the latent state
parameters.
%matplotlib inline
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
# Get some data
nile = sm.datasets.nile.load_pandas().data['volume']
nile.index = pd.DatetimeIndex(start="1871",end="1971", freq='y')
# Estimate the parameters of the model
mod = sm.tsa.UnobservedComponents(nile, 'llevel')
res = mod.fit()
# Get the simulation smoother and simulate
sim = res.model.simulation_smoother()
sim.simulate()
# Plot the filtered state and the state produced by the simulation smoother
pd.Series(sim.generated_obs[0,:]).plot(observed=False)
res.plot_components()
I would have expected state produced by the simulation smoother to look
similar to the estimated state from the fitted model. However, as far as I
can see this is just simulating a local level process with the given
parameters.
I have built statsmodels from the version on github (0.10.0 according to
the output of "conda list")
Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I am having trouble getting simulation smoothing of state space models
to work.
Suppose I have an UnobservedComponents model.
fitted = mod.fit()
sim = mod.simulation_smoother()
sim.simulate(initial_state_variates=fitted.smoothed_state[:,0],
pretransformed_variates=True)
However the resulting output does not seem to be conditional on the
observations. Is there a kwarg I need to pass to the simulation_smoother to
get this to work.
I have tried to figure this out from studying the source but to no avail.
Best wishes,
Alastair

There are no kwargs that you should have to pass.
It's hard to say more based on your e-mail, though. Can you describe
what you're trying to do in more detail, e.g. what you're expecting and
what's actually happening, and the version of Statsmodels that you're using?

Hi Alastair,
Thanks for posting that!
I think maybe you're not looking at the variable you want - I'm guessing
that you're looking for `sim.simulated_state`, which is the draw from the
joint posterior of the states.
Best,
Chad

[Chad Fulton]

Post by Alastair Heggie
Great, that works, not sure how I missed that...
For reference here's a notebook summarizing the use case I had in mind
https://gist.github.com/AliHeggie/3f2f56963a0e6893043d4116dfa3caed
Best,
Alastair

This is a very interesting notebook. I haven't used the NUTS sampler
before, and I notice that it provides a quite a bit different posterior for
the "obs_cov" than that obtained from Gibbs sampling, as computed in
http://www.chadfulton.com/topics/local_level_nile.html, see e.g. the trace
plot for \sigma_\varepsilon^2 (top left subplot in the graph right above
the "Expanded model: UnobservedComponents" heading).

(In fact, the posterior from an MH sampler with RW proposal looks quite a
bit like that from the Gibbs sampling too).

Best,
Chad

[Alastair Heggie]

I just did a few experiments to see what might be causing the difference
but now I'm more confused

updated gist:
https://gist.github.com/AliHeggie/3f2f56963a0e6893043d4116dfa3caed

My original notebook used a Gamma prior with alpha=2 and beta = 2/sigma
where sigma is the MLE estimate of the parameter.

Since the mean of a gamma distribution is alpha/beta the mean of this prior
is the MLE estimate of the variance I think. I don't know if it is poor
practice to use MLE estimates to inform your priors.

Anyway, one side effect this might be that the variance of the prior is too
low. So I tried using this method with alpha= 10 and the results seem
closer to your gibbs sampling notebook although not the same.

I then tried using the NUTS sampler with what I thought were your priors:
gamma(3,3). The results are very different with the mean of the posterior
for both parameters several orders of magnitude smaller. Although I now
notice you were using a inverse gamma in that notebook, so this might not
be very significant.

I'm not sure how to use a gibbs sampler with Pymc3 since
mc.step_methods.gibbs is a module.

Metropolis and HamiltonianMC seem to give similar results to NUTS but the
chains are less well mixed.

Post by Chad Fulton

Post by Alastair Heggie
Great, that works, not sure how I missed that...
For reference here's a notebook summarizing the use case I had in mind
https://gist.github.com/AliHeggie/3f2f56963a0e6893043d4116dfa3caed
Best,
Alastair

This is a very interesting notebook. I haven't used the NUTS sampler
before, and I notice that it provides a quite a bit different posterior for
the "obs_cov" than that obtained from Gibbs sampling, as computed in
http://www.chadfulton.com/topics/local_level_nile.html, see e.g. the
trace plot for \sigma_\varepsilon^2 (top left subplot in the graph right
above the "Expanded model: UnobservedComponents" heading).
(In fact, the posterior from an MH sampler with RW proposal looks quite a
bit like that from the Gibbs sampling too).
Best,
Chad