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!

No problem. I think the policy is to provide functions to the user that are defined on the ``MarkovChain`` type. That way if we introduce different kinds of Markov chains that require different functions to perform the same calculations, we can just overload the existing functions. E.g. ``mc_compute_stationary(dmc::DiscreteMarkovChain)``, ``mc_compute_statinary(cmc::ContinuousMarkovChain)``.

If we're going with GTH as the method in ``mc_compute_stationary`` then on an irreducible Markov chain, this will call ``gth_solve``. I kept ``lu_solve`` and ``eigen_solve`` hanging around in case we wanted to add an option to ``mc_compute_stationary`` that allows the user to pick the solution algorithm. 

Basically for this reason I think it makes sense to keep these subroutines hidden "behind the scenes".