MarkovChain simulate: Match with python version

[oyamad]

This is to propose matching the simulate method with that in the Python version (see QuantEcon/QuantEcon.py#146 (comment)). It works as follows:

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> simulate(mc, 5)
5-element Array{Int64,1}:
 2
 2
 1
 1
 2

julia> simulate(mc, 5, [1, 2])
5x2 Array{Int64,2}:
 1  2
 1  1
 2  1
 2  1
 2  2

julia> simulate(mc, 5, num_reps=3)
5x3 Array{Int64,2}:
 1  2  2
 2  1  1
 2  2  1
 1  2  2
 1  1  1

julia> simulate(mc, 5, 1, num_reps=3)
5x3 Array{Int64,2}:
 1  1  1
 1  1  2
 2  2  2
 2  2  2
 1  2  2

julia> simulate(mc, 5, [1, 2], num_reps=3)
5x6 Array{Int64,2}:
 1  2  1  2  1  2
 2  2  1  1  2  1
 2  1  2  2  1  2
 2  2  1  2  1  2
 2  2  2  2  2  2

As discussed in #52 (comment), it uses searchsortedfirst (via DiscreteRV). The method using searchsortedfirst (simulate) seems faster than the previous one using Categorical (mc_sample_path) for medium- and large-size Markov chains; see mc_simulate.jl.ipynb.

In passing, Numba seems to do a good job; see mc_simulate.py.ipynb.

I am not sure if I properly exploited the multiple dispatch system of Julia. Any suggestion/correction etc will be welcome!

[sglyon]

No need for the Union and nothing here. Just set a default argument of num_reps=1

[oyamad]

Done.

[oyamad]

Reverted.

[sglyon]

Looks great overall, thanks for adding this.

I have another request: can we document how the columns are ordered when num_reps > 1?

[oyamad]

"If init is a vector, the sample path from the j-th repetition of the simulation with initial state init[i] is stored in the i*j-th (j-1)*num_reps+i-th column of the matrix X."

Is this what you mean?

[sglyon]

Sorry this has sat for so long.

@oyamad what's the status on this, do you recall?

[oyamad]

@spencerlyon2 Sorry I had forgotten about this PR.

I now remembered why I did num_reps::Union{Int, Void}=nothing in the first simulate method; it is called from the other methods with num_reps=nothing. There may be a good way to avoid it, but I left it as is. Anyway, the issue is whether it is a good implementation that if num_reps is not supplied, the output is a 1-dimensitonal array, while if it is supplied, the output is 2-dimensional even when num_reps=1.

[sglyon]

Ahh that's right thanks @oymad.

The behavior you described will be type stable, so I think it is a good choice. I'm happy to merge as is.

Thanks

[oyamad]

@spencerlyon2 Thanks for taking care of the tests.

[sglyon]

You're welcome, thanks again for the quality work!