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!

I'm not sure that these type procedures and dimensions checking are strictly necessary since on operations between two `Int`'s the result is automatically promoted to a float, e.g. ``1/2 = 0.5``, and the `MarkovChain` type already ensures correct dimensions. Does the algorithm run okay without these portions?