jstac/quant-econ — pull request #79 by oyamad (2014-09-10 05:29:14+00:00)
Title: MARKOV: More efficient implementation for computing stationary distributions

Hi,
related to #18, I wrote some code that implements an algorithm specialized in computing stationary distributions (instead of using a general routine to solve the eigenvalue problem), called the Grassmann-Taksar-Heyman (GTH) algorithm (see for example [these slides](http://www.sti.uniurb.it/events/sfm07pe/slides/Stewart_1.pdf), pages 46-81). It is more straightforward than the current procedure, and it turned out that it is accurate enough to give (almost) correct answers even for KMR matrices (without mpmath).

stochmatrix.py:
* `gth_solve` solves for a nonzero solution to x (P - I) = 0 by the GTH algorithm for an *irreducible* stochastic matrix P.
* For a possibly reducible matrix P, `StochMatrix` decomposes P into irreducible components (where [`scipy.sparse.csgraph.connected_components`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csgraph.connected_components.html) is used to find communication classes).
* `stationary_dists` passes these irreducible components to `gth_solve` to obtain all the stationary distributions, one for each irreducible component. (They are contained in the output array as its *rows*.)

mc_tools.py:
* `mc_compute_stationary` is replaced with the above `stationary_dists`.
* Attributes/methods of `StochMatrix` are made available in `DMarkov`, such as `comm_classes` and `rec_classes`.
* In passing, I changed the names of methods of `DMarkov` from "mc_compute_stationary" and "mc_sample_path" to "compute_stationary" and "simulate".

Please see [stochmatrix_demo01.ipynb](http://nbviewer.ipython.org/github/oyamad/stochmatrix/blob/master/stochmatrix_demo01.ipynb) and [dmarkov_demo01.ipynb](http://nbviewer.ipython.org/github/oyamad/stochmatrix/blob/master/dmarkov_demo01.ipynb) for some more details.

I also prepared test_stochmatrix.py (and edited test_mc_tools.py).

There are some issues/questions:

1. I am not very experienced in python (and in programming in general). In particular, in `StochMatrix` I use decorators, but I am not sure if they are properly used, as well as if my choices between attribute and method in implementing each functionality are appropriate.

2. `stationary_dists` is slower than the current `mc_compute_stationary`, but I hope it is acceptable.

3. In `DMarkov`:
1. The argument `pi_0` is not used. Is it for any future use?
2. It has an attribute `stationary_dists` in addition to the method to compute stationary distributions. Is it only for the `__repr__` method?
3. It checks if the rows of the input matrix sum exactly to one. I am afraid it is too strict, and it might be better to allow some tolerance.

4. Test is missing for `mc_sample_path`. I could add one.