Refactor markov

src/QuantEcon.jl

@@ -41,8 +41,11 @@ export
 
 # mc_tools
     MarkovChain,
-    mc_compute_stationary, mc_sample_path, mc_sample_path!,
-    period, is_irreducible, is_aperiodic, recurrent_classes, communication_classes, simulate,
+    mc_compute_stationary,
+    simulate, simulate!, simulate_values, simulate_values!,
+    simulation, value_simulation,
+    period, is_irreducible, is_aperiodic, recurrent_classes,
+    communication_classes, n_states,
 
 # gth_solve
     gth_solve,
@@ -115,9 +118,10 @@ include("util.jl")
 ##### includes
 include("arma.jl")
 include("compute_fp.jl")
-include("markov_approx.jl")
-include("mc_tools.jl")
+include("markov/markov_approx.jl")
+include("markov/mc_tools.jl")
 include("markov/ddp.jl")
+include("markov/random_mc.jl")
 include("discrete_rv.jl")
 include("ecdf.jl")
 include("estspec.jl")
@@ -130,6 +134,5 @@ include("matrix_eqn.jl")
 include("robustlq.jl")
 include("quad.jl")
 include("quadsums.jl")
-include("random_mc.jl")
 
 end # module

src/discrete_rv.jl

@@ -25,17 +25,22 @@ vector of probabilities given by q.
 
 ##### Fields
 
-- `q::Vector{T<:Real}`: A vector of non-negative probabilities that sum to 1
-- `Q::Vector{T<:Real}`: The cumulative sum of q
+- `q::AbstractVector`: A vector of non-negative probabilities that sum to 1
+- `Q::AbstractVector`: The cumulative sum of q
 """
-type DiscreteRV{T<:Real}
-    q::Vector{T}
-    Q::Vector{T}
-    DiscreteRV(x::Vector{T}) = new(x, cumsum(x))
+type DiscreteRV{TV1<:AbstractVector, TV2<:AbstractVector}
+    q::TV1
+    Q::TV2
+    function DiscreteRV(q, Q)
+        abs(Q[end] - 1) > 1e-10 && error("q should sum to 1")
+        new(q, Q)
+    end
 end
 
-# outer constructor so people don't have to put {T} themselves
-DiscreteRV{T<:Real}(x::Vector{T}) = DiscreteRV{T}(x)
+function DiscreteRV{TV<:AbstractVector}(q::TV)
+    Q = cumsum(q)
+    DiscreteRV{TV,typeof(Q)}(q, Q)
+end
 
 """
 Make a single draw from the discrete distribution
@@ -63,4 +68,4 @@ Make multiple draws from the discrete distribution represented by a
 
 - `out::Vector{Int}`: `k` draws from `d`
 """
-draw{T}(d::DiscreteRV{T}, k::Int) = Int[draw(d) for i=1:k]
+draw(d::DiscreteRV, k::Int) = Int[draw(d) for i=1:k]

src/markov_approx.jl → src/markov/markov_approx.jl

@@ -34,8 +34,8 @@ should span
 
 ##### Returns
 
-- `y::Vector{Real}` : Nodes in the state space
-- `Π::Matrix{Real}` Matrix transition probabilities for Markov Process
+- `mc::MarkovChain{Float64}` : Markov chain holding the state values and
+transition matrix
 
 """
 function tauchen(N::Integer, ρ::Real, σ::Real, μ::Real=0.0, n_std::Integer=3)
@@ -75,7 +75,11 @@ function tauchen(N::Integer, ρ::Real, σ::Real, μ::Real=0.0, n_std::Integer=3)
 
     y .+= μ / (1 - ρ) # center process around its mean (wbar / (1 - rho))
 
-    return y, Π
+    # renormalize. In some test cases the rows sum to something that is 2e-15
+    # away from 1.0, which caused problems in the MarkovChain constructor
+    Π = Π./sum(Π, 2)
+
+    MarkovChain(Π, y)
 end
 
 
@@ -96,8 +100,8 @@ where ε_t ~ N (0, σ^2)
 
 ##### Returns
 
-- `y::Vector{Real}` : Nodes in the state space
-- `Θ::Matrix{Real}` Matrix transition probabilities for Markov Process
+- `mc::MarkovChain{Float64}` : Markov chain holding the state values and
+transition matrix
 
 """
 function rouwenhorst(N::Integer, ρ::Real, σ::Real, μ::Real=0.0)
@@ -107,14 +111,15 @@ function rouwenhorst(N::Integer, ρ::Real, σ::Real, μ::Real=0.0)
     ψ = sqrt(N-1) * σ_y
     m = μ / (1 - ρ)
 
-    return rouwenhorst(p, p, m, ψ, N)
+    state_values, p = _rouwenhorst(p, p, m, ψ, N)
+    MarkovChain(p, state_values)
 end
 
-function rouwenhorst(p::Real, q::Real, m::Real, Δ::Real, n::Integer)
+function _rouwenhorst(p::Real, q::Real, m::Real, Δ::Real, n::Integer)
     if n == 2
-        return Real[m-Δ, m+Δ], [p 1-p; 1-q q]
+        return [m-Δ, m+Δ],  [p 1-p; 1-q q]
     else
-        _, θ_nm1 = rouwenhorst(p, q, m, Δ, n-1)
+        _, θ_nm1 = _rouwenhorst(p, q, m, Δ, n-1)
         θN = p    *[θ_nm1 zeros(n-1, 1); zeros(1, n)] +
              (1-p)*[zeros(n-1, 1) θ_nm1; zeros(1, n)] +
              q    *[zeros(1, n); zeros(n-1, 1) θ_nm1] +