How do you do interventions on a statespace model?

I’ve calibrated a statespace model and want to see how its behaviour changes when I modify some of the model parameters. I thought I might be able to use the do function, but if that’s possible I haven’t been able to figure out the right syntax to use.

CC @jessegrabowski, this may require some understanding of the internals of Statespace models (if you’re using pymc-extras) or there may be some built-in functionality (heard the shock term being thrown around before)

Ideally, you would be able to use the do operator. Currently I think it’s not possible, and you instead have to jump through a few hoops. I’m going to give you the full rundown of how statespace works, then at the end tell you what I’d do in your case. My feelings won’t be hurt if you skip to the end.

  • The statespace model defines a list of parameters that are expected as random variables, plus a construction method for where all those parameters go in the various statespace matrices
  • When you call build_statespace_graph inside a pymc model context, three internal helpers are called: ss_mod._insert_random_variables(), ss_mod._save_exogenous_data_info() ss_mod._insert_data_variables(). The second two are just bookkeeping helpers so that you can do pm.set_data on regression-type problems, let’s set them aside.
  • _insert_random_variables takes the RVs found in the current PyMC model context and matches them by name with expected parameter names. For example, if the statespace model asks for initial_level_trend, it looks for an RV of that name and plugs it in everywhere it is supposed to go.
  • Ok “plugs into” what? The statespace model maintains two sets of statespace matrices. self.ssm maintains the “raw” StatespaceRepresentation class. This holds the parameter → matrix mappings, and never gets mutated. self.subbed_ssm is what gets mutated when you call _insert_random_variables()
  • There is a helper function get_statespace_matrices that returns the substituted x0, P0, c, d, T, Z, R, H, Q objects.

Now, all this matching only cares about name, it does not care about where a thing comes from. So you can, for example, do pm.Deterministic('initial_level_trend', [3., 0.1], dims=whatever) and it will insert the constant vector [3., 0.1] in that slot. It doesn’t care.

Now, in your case that doesn’t help because you already sampled and you have an idata. You want to call a post-estimation function like forecast, impulse_response_function, or sample_conditional_posterior. These functions work by internally setting up a dummy PyMC model with fake random variables, then calling self._insert_random_variables() to plug them in, then sampling from that dummy model using the draws from idata. This is the normal pm.sample_posterior_predictive machine explained e.g. here. Technically this is the correct moment to do pm.do. Ideally, you want to call it on this dummy model, then PyMC should simply do the right thing. Statespace might not though – if pm.do removes the named random variables from the model, the statespace model name-matching will be confused and throw an error complaining than an expected parameter isn’t there.

So what can you do? One option would be to add a a replacements argument to the post-estimation methods that would let you map parameters to fixed values. It would do the trick but it would make the signatures of these methods even more horrible.

Perhaps a better option would be to let you call these post estimation methods inside a pymc model context, in which case we take the open model context as the “dummy model” and just set up the forecast/irf/filtering problem taking that model as given. Maybe better, but would be a bit of a “hidden feature”.

If you’re not interested in boiling the ocean and just want to do some experiments, I would say the easiest way forward for you right now is to copy the idata and replace the draws of the variables you want to intervene on with the values you want to test. When you call the post-estimation methods, they take the idata as given, so whatever you do there will be correctly propagated. It’s not pretty but it’s the best i can offer.