Merge pull request #26 from ZacCranko/add-gth-support

Add gth support to mc_compute_stationary

src/QuantEcon.jl

@@ -110,8 +110,7 @@ include("discrete_rv.jl")
 include("distributions.jl")
 include("ecdf.jl")
 include("estspec.jl")
-include("gth_solve.jl")
-include("kalman.jl")
+include("mc_tools.jl")
 include("kalman.jl")
 include("lae.jl")
 include("lqcontrol.jl")

src/mc_tools.jl

@@ -41,20 +41,23 @@ function Base.show(io::IO, mc::MarkovChain)
 end
 
 # function to solve x(P-I)=0 by eigendecomposition
-function eigen_solve(p::Matrix)
-    ef = eigfact(p)
+function eigen_solve{T}(p::Matrix{T})
+    ef = eigfact(p')
     isunit = map(x->isapprox(x,1), ef.values)
-    x = ef.vectors[:, isunit]
-    x ./= norm(x,1) # normalisation
-    any(x .< 0) && throw("something has gone wrong with the eigen solve")
-    abs(x) # avoid entries like '-0' appearing
+    x = real(ef.vectors[:, isunit])
+    x ./= sum(x,1) # normalisation
+    for i = 1:length(x)
+        x[i] = isapprox(x[i],zero(T)) ? zero(T) : x[i]
+    end
+    any(x .< 0) && warn("something has gone wrong with the eigen solve")
+    x
 end
 
 # function to solve x(P-I)=0 by lu decomposition
 function lu_solve{T}(p::Matrix{T})
     n,m = size(p)
     x   = vcat(Array(T,n-1),one(T))
-    u   = lufact(p - one(p))[:U]
+    u   = lufact(p' - one(p))[:U]
     for i = n-1:-1:1 # backsubstitution
         x[i] = -sum([x[j]*u[i,j] for j=i:n])/u[i,i]
     end
@@ -66,6 +69,40 @@ function lu_solve{T}(p::Matrix{T})
     x
 end
 
+gth_solve{T<:Integer}(A::Matrix{T}) = gth_solve(float64(A))
+
+function gth_solve{T<:Real}(A::AbstractMatrix{T})
+    A1 = copy(A)
+    n = size(A1, 1)
+    x = zeros(T, n)
+
+    # === Reduction === #
+    for k in 1:n-1
+        scale = sum(A1[k, k+1:n])
+        if scale <= 0
+            # There is one (and only one) recurrent class contained in
+            # {1, ..., k};
+            # compute the solution associated with that recurrent class.
+            n = k
+            break
+        end
+        A1[k+1:n, k] ./= scale
+
+        for j in k+1:n, i in k+1:n
+            A1[i, j] += A1[i, k] * A1[k, j]
+        end
+    end
+
+    # === Backward substitution === #
+    x[end] = 1
+    for k in n-1:-1:1, i in k+1:n
+        x[k] += x[i] * A1[i, k]
+    end
+
+    # === Normalization === #
+    x / sum(x)
+end
+
 # find the reducible subsets of a markov chain
 function irreducible_subsets(mc::MarkovChain)
     p = bool(mc.p)
@@ -99,20 +136,23 @@ end
 # calculate the stationary distributions associated with a N-state markov chain
 # output is a N x M matrix where each column is a stationary distribution
 # currently using lu decomposition to solve p(P-I)=0
-function mc_compute_stationary(mc::MarkovChain)
+function mc_compute_stationary(mc::MarkovChain; method=:gth)
+    solvers = Dict([:gth => gth_solve, :lu => lu_solve, :eigen => eigen_solve])
+    solve = solvers[method]
+
     p,T = mc.p,eltype(mc.p)
     classes = irreducible_subsets(mc)
 
     # irreducible mc
-    length(classes) == 1 && return lu_solve(p')
+    length(classes) == 1 && return solve(p)
 
     # reducible mc
     stationary_dists = Array(T,n_states(mc),length(classes))
     for i = 1:length(classes)
         class  = classes[i]
         dist   = zeros(T,n_states(mc))
         temp_p = p[class,class]
-        dist[class] = lu_solve(temp_p')
+        dist[class] = solve(temp_p)
         stationary_dists[:,i] = dist
     end
     return stationary_dists
@@ -148,5 +188,7 @@ end
 # simulate markov chain starting from some initial value. In other words
 # out[1] is already defined as the user wants it
 function mc_sample_path!(mc::MarkovChain, samples::Vector)
-    samples = mc_sample_path(mc,samples[1],samples)
+    length(samples) < 2 &&
+        throw(ArgumentError("samples vector must have length greater than 2"))
+    samples = mc_sample_path(mc,samples[1],length(samples)-1)
 end

test/runtests.jl

@@ -7,7 +7,6 @@ include("test_compute_fp.jl")
 include("test_discrete_rv.jl")
 include("test_ecdf.jl")
 include("test_estspec.jl")
-include("test_gth_solve.jl")
 include("test_kalman.jl")
 include("test_lae.jl")
 include("test_lqcontrol.jl")

test/test_mc_tools.jl

@@ -20,13 +20,22 @@ P5 = [
      ]
 P5_stationary = hcat([1/2, 1/2, 0, 0, 0, 0],[0, 0, 0, 1/3, 1/3, 1/3])
 
+P6 = [2//3 1//3; 1//4 3//4]  # Rational elements
+P6_stationary = [3//7, 4//7]
+
+P7 = [1 0; 0 1]
+P7_stationary = [1 0;0 1]
+
 d1 = MarkovChain(P1)
 d2 = MarkovChain(P2)
 d3 = MarkovChain(P3)
 d4 = MarkovChain(P4)
 d5 = MarkovChain(P5)
+d6 = MarkovChain(P6)
+d7 = MarkovChain(P7)
 
-function KMR_Markov_matrix_sequential(N, p, epsilon)
+
+function kmr_markov_matrix_sequential(N, p, epsilon)
     """
     Generate the Markov matrix for the KMR model with *sequential* move
 
@@ -52,14 +61,40 @@ function KMR_Markov_matrix_sequential(N, p, epsilon)
     return P
 end
 
-facts("Testing mc_tools.jl") do
+P8 = kmr_markov_matrix_sequential(27, 1/3, 1e-2)
+P9 = kmr_markov_matrix_sequential(3, 1/3, 1e-14)
 
+d8 = MarkovChain(P8)
+d9 = MarkovChain(P9)
+
+tol = 1e-15
+
+facts("Testing mc_tools.jl") do
     context("test mc_compute_stationary using exact solutions") do
         @fact mc_compute_stationary(d1) => eye(3)[:, [1, 3]]
         @fact mc_compute_stationary(d2) => roughly([0, 9/14, 5/14])
         @fact mc_compute_stationary(d3) => roughly([1/4, 3/4])
         @fact mc_compute_stationary(d4) => eye(2)
         @fact mc_compute_stationary(d5) => roughly(P5_stationary)
+        @fact mc_compute_stationary(d6) => P6_stationary
+        @fact mc_compute_stationary(d7) => roughly(P7_stationary)
+    end
+
+    context("test gth_solve with KMR matrices") do
+        for d in [d8,d9]
+            x = mc_compute_stationary(d)
+
+            # Check elements sum to one
+            @fact sum(x) => roughly(1; atol=tol)
+
+            # Check elements are nonnegative
+            for i in 1:length(x)
+                @fact x[i] => greater_than_or_equal(-tol)
+            end
+
+            # Check x is a left eigenvector of P
+            @fact vec(x'*d.p) => roughly(x; atol=tol)
+        end
     end
 
     context("test MarkovChain throws errors") do
@@ -69,7 +104,4 @@ facts("Testing mc_tools.jl") do
     end
 end  # facts
 
-end  # module
-
-# TODO: P = KMR_Markov_matrix_sequential(27, 1/3, 1e-2) will fail without
-#       arbitrary precision linear algebra. Need to wait on Juila for this
+end  # module