ChadFulton
commented
7 years ago Hi Mathieu,
First, just to make sure, the pymar library should be considered as obsolete - a much improved version is available in Statsmodels v0.8 - v0.8 is not released yet, but there is a release candidate on pip. Alternatively Statsmodels master can be used.
The equation you're referencing is a state space model. In that paper it is estimated via Bayesian methods, and the Gibbs sampling step essentially requires support for a Markov switching vector autoregression.
The current Statsmodels switching models are all univariate, but there should not be any difficulty expanding to a VAR (the Hamilton filter / Kim smoother only requires the computed likelihoods). This is in my todo list (no timeline though, it's whenever I get back to that project); the real impediment isn't writing the switching VAR code, it's that getting unit tests together always takes some time.
To implement it, I suggest extending the MarkovSwitching class (see the MarkovRegression and MarkovAutoregression classes for examples) with e.g. a MarkovVAR class. The only thing you really need to implement is the _conditional_likelihoods method, which computes the likelihood for each observation conditional on the parameters and the current and, since it's a VAR, past regimes. See MarkovAutoregression for the univariate case. For example, for a VAR(2), you need to condition on the current regime as well as the past two regimes. So the _conditional_likelihoods method will end up returning a 2x2x2xnobs array, where the (i,j,k,s) element is the likelihood of the s-th observation when s is in regime i, s-1 is in regime j, and s-2 is in regime s-3.
We're not quite set up to perform Gibbs sampling in a state space framework; there is a simulation smoother available in a pull request, that will be merged after v0.8 is released (of course you could use the PR's branch right now).
If you want to estimate the state space model by maximum likelihood, then it's a different procedure. You need to use Kim's approximate version of the Hamilton filter; that would more difficult to implement. There was a GSoC project on that subject this last summer, but it was in Python and its performance is not so good yet. I have improved versions in Cython, but they're not ready yet.
How hard to you think it would be to generalize your model to one where you have multiple time series all sharing the same switching factor (like Equation (1) in http://pages.uoregon.edu/jpiger/research/published-papers/chauvet-and-piger_2008_jour.pdf )
And do you have any tips to someone who might give this a try?