QuantEcon/QuantEcon.jl — pull request #23 by oyamad (2014-11-24 07:31:16+00:00)
Title: Implement the GTH algorithm

This is to import [gth_solve.py](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/gth_solve.py) to the Julia part (I also imported the tests as well).
(I have had a brief discussion with @ZacCranko on this matter.)

1.
This is why I think `gth_solve` is necessary also for the Julia side:

I guess that the source of the problem is the same as in NumPy and SciPy, but it seems the routine (probably `lufact`) from the Julia standard library loses precision for some inputs with very small float numbers.

Here's an example:

```julia
for i in 12:18
a = 1e-1 ^ i
P = [1-a a; a/2 1-a/2]
mc = MarkovChain(P)
print("a = 1e-", i, "\n")
print(mc_compute_stationary(mc), "\n")
end
```

The output is:

```julia
a = 1e-12
[0.3333530884328846,0.6666469115671155]
a = 1e-13
[0.33326424872384786,0.6667357512761523]
a = 1e-14
[0.33315566870167473,0.6668443312983253]
a = 1e-15
[0.3569592970777264,0.6430407029222737]
a = 1e-16
[0.3105159921038858,0.6894840078961141]
a = 1e-17
[0.0,1.0]
a = 1e-18
[0.0,1.0]
```

The correct answer is [1/3, 2/3] for any value of `a`.

A problem occurs also for "KMR Markov matrices" with small epsilon.

The "GTH" algorithm is robust to this kind of inputs, as it does not involve subtraction. This is why I included "gth_solve.py" into the python side, and I believe this is also necessary for the julia side.

2.
I de-vectorized the python version following the general advice for Julia. Surprisingly, this de-vectorized version is 10x faster than the direct translation of the python version using array indexing for some random 1000x1000 matrices.

3.
Docstring is yet to be imported.

4.
I tried to follow a suggestion by @ZacCranko and manipulated the code, and then it appears the routine now returns an output with rational number elements for a matrix input with rational number elements:

```julia
P = [2//3 1//3; 1//4 3//4]
gth_solve(P)
```

Output:

```julia
2-element Array{Rational{Int64},1}:
3//7
4//7
```

But I don't think I really understand what is going on..

Anyway, this is effectively my first code in Julia, and definitely needs reviews by experienced Julia programmers... Thanks!