New features for mc_tools

[ZacCranko]

• Moved to LightGraphs for graph calculations
• New functions added

recurrent_classes(mc::MarkovChain)
communication_classes(mc::MarkovChain)
is_irreducible(mc::MarkovChain)
period(mc::MarkovChain)

Tests for these new functions are not included in this pull. I'm going to change up the tests file, so I'm submitting this pull and making those modifications separately, but all the old tests pass fine with the new code.

[oyamad]

@ZacCranko Thanks, it looks very nice.

1. It will be helpful to have is_aperiodic as well:

   is_aperiodic(mc::MarkovChain) = period(DiGraph(mc.p)) == 1
   

   or add is_aperiodic to LightGraphs.jl and do

   is_aperiodic(mc::MarkovChain) = is_aperiodic(DiGraph(mc.p))
   

2. In mc_tools.py, the period is also defined for a reducible Markov chain, to be the LCM of the periods of the recurrent classes. I think it is better to adopt the same definition between the two versions.

3. Also remaining is cyclic_classes.

4. For method in mc_compute_stationary, by what criterion should the user choose among gth, lu, and eigen (what is the point of having three methods)?

[sglyon]

We should restore the inner+outer constructor pair. The reason is that we never want to be able to create an instance of MarkovChain without these 3 conditions holding -- these are the invariants referred to in the inner constructor section of the manual

[oyamad]

WIth Version 0.3.10, I got the following error:

julia> using QuantEcon

julia> mc = MarkovChain([0.4 0.6; 0.2 0.8])
Discrete Markov Chain
stochastic matrix:
2x2 Array{Float64,2}:
 0.4  0.6
 0.2  0.8


julia> mc_compute_stationary(mc)
ERROR: error compiling mc_compute_stationary: error compiling __mc_compute_stationary#140__: unsupported or misplaced expression => in function __mc_compute_stationary#140__

It works fine with Version 0.4.0.

[sglyon]

Thanks for the note. We'll take care of it

[sglyon]

This line needs @compat to work on 0.3 -- as @oyamad noted

[sglyon]

We should also make sure that the method argument is typed as ::Symbol (anything else will throw an error when we try to access the dict)

(we should restore exactly what was here originally -- it reverts both changes)

[ZacCranko]

Thanks guys, I believe I've addressed everything except for the import question, and cyclic_components which I'm not gonna tackle just yet

[oyamad]

@ZacCranko Thanks!

• mc_compute_stationary returns a 1-dimensional array if mc has a unique recurrent class (in which case there is a unique stationary state) and a 2-dimensional array otherwise, while it is always 1-dimensional in the Python version. Would there be any problem having different behaviors?
• What about the methods for mc_compute_stationary? What is your intention?
• This is not about this PR itself, but I think mc_sample_path(mc, init, sample_size) should return an array of length sample_size, not sample_size+1.

[ZacCranko]

I don't really have much of an opinion on the solution methods. Though since this is a pedagogical tool, there is probably utility in allowing the the user to compare GTH and the other methods, especially if we, in the future, were to have a lecture on approaches to compute the stationary distribution of a Markov chain. This would be consistent with the fact that we're exporting gth_solve. Any thoughts @spencerlyon2?

[oyamad]

I see, thanks. Makes sense for education purposes.

[ZacCranko]

Thanks for that, good to know! Comments fixed.

[oyamad]

Another issue: Why does mc_compute_stationary return stationary distributions as columns, whereas MarkovChain accepts a row-stochastic matrix, for which a stationary distribution is a left-eigenvector? I find this somewhat inconsistent.

[ZacCranko]

I'm happy to have mc_compute_stationary return what ever anyone would like. However since Julia is column-major it would be more efficient computationally (for accessing) if vectors are columns.

[oyamad]

Right, then we might consider having MarkovChain accept a column-stochastic matrix. Mathematically, there is no reason for one of row-stochastic or column-stochastic to be better than the other.

[ZacCranko]

The only convenience with the current orientation of the stochastic matrix is that it coincides with the definition of an adjacency matrix.

[oyamad]

Ah that's a good point.

Would it be costly to transpose the column-stochastic input and pass it to DiGraph?

[ZacCranko]

Computing the transpose of a matrix is O(n), but this is somewhat unrelated to what this PR is addressing. If you really feel strongly about this I think it would be best to open an issue.

[oyamad]

Yes, that's right.

[jstac]

Sorry if I get something wrong here as I'm just going by the thread, but I would recommend that we return the stationary distributions -- and all other distributions -- as row vectors, either by transposing them as they are passed out or some other approach. It's entirely conventional in the field of Markov chains that distributions are row vectors. Similarly, the matrix should be such that P[i, j] means prob of moving from state i to state j. (In fact that's why distributions end up being row vectors that we multiply from the left.) So the challenge is to implement this standard interface in the most efficient way.

[sglyon]

What if we return the stationary distributions as a Vector{Vector{Float64} (Or whatever T<:Real we computed them as instead of Float64)? Then we don't have to worry about orientation

[jstac]

I guess a flat array is fine.

[ZacCranko]

Ahoy hoy, is this ready to be merged?

[sglyon]

Hey sorry, I kinda forgot we hadn't merged this.

I would love to see some tests of the new methods. i realize you want to change the tests. Will that happen soon?

[ZacCranko]

Yep, that'll happen pretty quick.

[sglyon]

Great, I'll break my rule of merging without tests this time...

Thanks for the work on this.