From cacb9ef28aed17c22beab81f95119aadd56b7f1e Mon Sep 17 00:00:00 2001 From: Chase Date: Tue, 23 Feb 2016 11:47:59 -0500 Subject: [PATCH 1/9] Move markovapprox and mctools to markov folder. --- src/{ => markov}/markov_approx.jl | 0 src/{ => markov}/mc_tools.jl | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename src/{ => markov}/markov_approx.jl (100%) rename src/{ => markov}/mc_tools.jl (100%) diff --git a/src/markov_approx.jl b/src/markov/markov_approx.jl similarity index 100% rename from src/markov_approx.jl rename to src/markov/markov_approx.jl diff --git a/src/mc_tools.jl b/src/markov/mc_tools.jl similarity index 100% rename from src/mc_tools.jl rename to src/markov/mc_tools.jl From 52f8b8fbda2c6f1b7750863d776f9899d3a6f3c0 Mon Sep 17 00:00:00 2001 From: Chase Date: Tue, 23 Feb 2016 12:33:59 -0500 Subject: [PATCH 2/9] Add folder structure for markov in QE.jl --- src/QuantEcon.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl index f8ffa92f..8c064560 100644 --- a/src/QuantEcon.jl +++ b/src/QuantEcon.jl @@ -115,8 +115,8 @@ 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("discrete_rv.jl") include("ecdf.jl") From 430ba61f1da335c4e90a67911a770dd9fe073f83 Mon Sep 17 00:00:00 2001 From: Chase Date: Tue, 23 Feb 2016 13:46:42 -0500 Subject: [PATCH 3/9] Abstract types in discrete rv. --- src/discrete_rv.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/discrete_rv.jl b/src/discrete_rv.jl index eea532d2..09a74c03 100644 --- a/src/discrete_rv.jl +++ b/src/discrete_rv.jl @@ -29,13 +29,13 @@ vector of probabilities given by q. - `Q::Vector{T<:Real}`: The cumulative sum of q """ type DiscreteRV{T<:Real} - q::Vector{T} - Q::Vector{T} - DiscreteRV(x::Vector{T}) = new(x, cumsum(x)) + q::AbstractVector{T} + Q::AbstractVector{T} + DiscreteRV(x::AbstractVector{T}) = new(x, cumsum(x)) end # outer constructor so people don't have to put {T} themselves -DiscreteRV{T<:Real}(x::Vector{T}) = DiscreteRV{T}(x) +DiscreteRV{T<:Real}(x::AbstractVector{T}) = DiscreteRV{T}(x) """ Make a single draw from the discrete distribution From ddd2baa537c0a174c28d2e03eaf81368dfcac84c Mon Sep 17 00:00:00 2001 From: Chase Date: Tue, 23 Feb 2016 13:47:14 -0500 Subject: [PATCH 4/9] Abstract types in mc_tools. --- src/markov/mc_tools.jl | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/markov/mc_tools.jl b/src/markov/mc_tools.jl index 8d5da5f5..82d44fb8 100644 --- a/src/markov/mc_tools.jl +++ b/src/markov/mc_tools.jl @@ -36,6 +36,7 @@ Steady State Distributions for Markov Chains, " Operations Research (1985), University Press, 2009. """ +@inline check_stochastic_matrix(P) = maxabs(sum(P, 2) - 1) < 1e-15 ? true : false """ Finite-state discrete-time Markov chain. @@ -49,18 +50,19 @@ transitions. - `p::Matrix` The transition matrix. Must be square, all elements must be positive, and all rows must sum to unity """ -type MarkovChain{T<:Real} - p::Matrix{T} # valid stochastic matrix +type MarkovChain{T, TM<:AbstractMatrix} + p::TM # valid stochastic matrix - - function MarkovChain(p::Matrix) + function MarkovChain(p::AbstractMatrix) n, m = size(p) - if n != m + if eltype(p) != T + throw(ArgumentError("Types must be consistent with given matrix")) + elseif n != m throw(ArgumentError("stochastic matrix must be square")) elseif any(p.<0) throw(ArgumentError("stochastic matrix must have nonnegative elements")) - elseif any(x->!isapprox(sum(x), 1), [p[i, :] for i = 1:n]) + elseif !check_stochastic_matrix(p) throw(ArgumentError("stochastic matrix rows must sum to 1")) end @@ -69,7 +71,7 @@ type MarkovChain{T<:Real} end # Provide constructor that infers T from eltype of matrix -MarkovChain{T<:Real}(p::Matrix{T}) = MarkovChain{T}(p) +MarkovChain(p::AbstractMatrix) = MarkovChain{eltype(p), typeof(p)}(p) "Number of states in the markov chain `mc`" n_states(mc::MarkovChain) = size(mc.p, 1) @@ -82,10 +84,10 @@ end # solve x(P-I)=0 by eigen decomposition "$(____solver_docs)" -function eigen_solve{T}(p::Matrix{T}) - ef = eigfact(p') - isunit = map(x->isapprox(x, 1), ef.values) - x = real(ef.vectors[:, isunit]) +function eigen_solve{T}(p::AbstractMatrix{T}) + vals, vecs = eigs(p') + isunit = map(x->isapprox(x, 1), vals) + x = real(vecs[:, isunit]) x ./= sum(x, 1) # normalisation for i = 1:length(x) x[i] = isapprox(x[i], zero(T)) ? zero(T) : x[i] @@ -96,7 +98,7 @@ end # solve x(P-I)=0 by lu decomposition "$(____solver_docs)" -function lu_solve{T}(p::Matrix{T}) +function lu_solve{T}(p::AbstractMatrix{T}) n, m = size(p) x = vcat(Array(T, n-1), one(T)) u = lufact(p' - one(p))[:U] @@ -112,10 +114,10 @@ function lu_solve{T}(p::Matrix{T}) end "$(____solver_docs)" -gth_solve{T<:Integer}(a::Matrix{T}) = gth_solve(convert(Array{Rational, 2}, a)) +gth_solve{T<:Integer}(a::AbstractMatrix{T}) = gth_solve(convert(AbstractArray{Rational, 2}, a)) # solve x(P-I)=0 by the GTH method -function gth_solve{T<:Real}(original::Matrix{T}) +function gth_solve{T<:Real}(original::AbstractMatrix{T}) a = copy(original) n = size(a, 1) x = zeros(T, n) From caff0fa9457162afe4673d1d83587b793949fffc Mon Sep 17 00:00:00 2001 From: Chase Date: Tue, 23 Feb 2016 13:47:49 -0500 Subject: [PATCH 5/9] Add room for sparse tests in mc_tools. --- test/test_mc_tools.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_mc_tools.jl b/test/test_mc_tools.jl index a6e64774..d94175a7 100644 --- a/test/test_mc_tools.jl +++ b/test/test_mc_tools.jl @@ -164,4 +164,8 @@ end @test vec(X[1, :]) == repmat(init, num_reps) end # testset + @testset "Sparse Matrix " begin + + end + end # testset From 800cc8ec9092c163518e127adbf7300274ace02f Mon Sep 17 00:00:00 2001 From: Spencer Lyon Date: Mon, 14 Mar 2016 09:48:42 -0400 Subject: [PATCH 6/9] final sparse MarkovChain support --- src/QuantEcon.jl | 5 +- src/markov/mc_tools.jl | 99 +++++++++++------------------------ src/{ => markov}/random_mc.jl | 0 test/test_mc_tools.jl | 68 ++++++++++++++++-------- 4 files changed, 81 insertions(+), 91 deletions(-) rename src/{ => markov}/random_mc.jl (100%) diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl index 8c064560..d22c99cd 100644 --- a/src/QuantEcon.jl +++ b/src/QuantEcon.jl @@ -42,7 +42,8 @@ export # mc_tools MarkovChain, mc_compute_stationary, mc_sample_path, mc_sample_path!, - period, is_irreducible, is_aperiodic, recurrent_classes, communication_classes, simulate, + period, is_irreducible, is_aperiodic, recurrent_classes, + communication_classes, simulate, n_states, # gth_solve gth_solve, @@ -118,6 +119,7 @@ include("compute_fp.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 +132,5 @@ include("matrix_eqn.jl") include("robustlq.jl") include("quad.jl") include("quadsums.jl") -include("random_mc.jl") end # module diff --git a/src/markov/mc_tools.jl b/src/markov/mc_tools.jl index 82d44fb8..b80a11ba 100644 --- a/src/markov/mc_tools.jl +++ b/src/markov/mc_tools.jl @@ -13,29 +13,6 @@ http://quant-econ.net/jl/finite_markov.html import LightGraphs: DiGraph, period, attracting_components, strongly_connected_components, is_strongly_connected -____solver_docs = """ -solve x(P-I)=0 using either an eigendecomposition, lu factorization, or an -algorithm presented by Grassmann-Taksar-Heyman (GTH) - -##### Arguments - -- `p::Matrix` : valid stochastic matrix - -##### Returns - -- `x::Matrix`: A matrix whose columns contain stationary vectors of `p` - -##### References - -The following references were consulted for the GTH algorithm - -- W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative Analysis and -Steady State Distributions for Markov Chains, " Operations Research (1985), -1107-1116. -- W. J. Stewart, Probability, Markov Chains, Queues, and Simulation, Princeton -University Press, 2009. - -""" @inline check_stochastic_matrix(P) = maxabs(sum(P, 2) - 1) < 1e-15 ? true : false """ @@ -76,45 +53,36 @@ MarkovChain(p::AbstractMatrix) = MarkovChain{eltype(p), typeof(p)}(p) "Number of states in the markov chain `mc`" n_states(mc::MarkovChain) = size(mc.p, 1) -function Base.show(io::IO, mc::MarkovChain) +function Base.show{T,TM}(io::IO, mc::MarkovChain{T,TM}) println(io, "Discrete Markov Chain") - println(io, "stochastic matrix:") + println(io, "stochastic matrix of type $TM:") println(io, mc.p) end -# solve x(P-I)=0 by eigen decomposition -"$(____solver_docs)" -function eigen_solve{T}(p::AbstractMatrix{T}) - vals, vecs = eigs(p') - isunit = map(x->isapprox(x, 1), vals) - x = real(vecs[:, 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 +""" +solve x(P-I)=0 using an algorithm presented by Grassmann-Taksar-Heyman (GTH) -# solve x(P-I)=0 by lu decomposition -"$(____solver_docs)" -function lu_solve{T}(p::AbstractMatrix{T}) - n, m = size(p) - x = vcat(Array(T, n-1), one(T)) - 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 - x ./= norm(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 lu solve") - x -end +##### Arguments + +- `p::Matrix` : valid stochastic matrix -"$(____solver_docs)" -gth_solve{T<:Integer}(a::AbstractMatrix{T}) = gth_solve(convert(AbstractArray{Rational, 2}, a)) +##### Returns + +- `x::Matrix`: A matrix whose columns contain stationary vectors of `p` + +##### References + +The following references were consulted for the GTH algorithm + +- W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative Analysis and +Steady State Distributions for Markov Chains, " Operations Research (1985), +1107-1116. +- W. J. Stewart, Probability, Markov Chains, Queues, and Simulation, Princeton +University Press, 2009. + +""" +gth_solve{T<:Integer}(a::AbstractMatrix{T}) = + gth_solve(convert(AbstractArray{Rational, 2}, a)) # solve x(P-I)=0 by the GTH method function gth_solve{T<:Real}(original::AbstractMatrix{T}) @@ -185,8 +153,6 @@ calculate the stationary distributions associated with a N-state markov chain ##### Arguments - `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix -- `;method::Symbol(:gth)`: One of `gth`, `lu`, and `eigen`; specifying which -of the three `_solve` methods to use. ##### Returns @@ -194,24 +160,22 @@ of the three `_solve` methods to use. distribution of `mc.p` """ -function mc_compute_stationary{T}(mc::MarkovChain{T}; method::Symbol=:gth) - solvers = Dict(:gth => gth_solve, :lu => lu_solve, :eigen => eigen_solve) - solve = solvers[method] - +function mc_compute_stationary{T}(mc::MarkovChain{T}) recurrent = recurrent_classes(mc) # unique recurrent class - length(recurrent) == 1 && return solve(mc.p) + length(recurrent) == 1 && return gth_solve(mc.p) # more than one recurrent classes stationary_dists = Array(T, n_states(mc), length(recurrent)) for (i, class) in enumerate(recurrent) dist = zeros(T, n_states(mc)) - dist[class] = solve(mc.p[class, class]) # recast to correct dimension + dist[class] = gth_solve(mc.p[class, class]) # recast to correct dimension stationary_dists[:, i] = dist end return stationary_dists end + """ Simulate a Markov chain starting from an initial state @@ -257,7 +221,7 @@ probability of being in seach state in the initial period """ function mc_sample_path(mc::MarkovChain, - init::Vector = vcat(1, zeros(n_states(mc)-1)), + init::Vector=vcat(1, zeros(n_states(mc)-1)), sample_size::Int=1000; burn::Int=0) # ensure floating point input for Categorical() init = convert(Vector{Float64}, init) @@ -294,7 +258,6 @@ function mc_sample_path!(mc::MarkovChain, samples::Array) for t=2:length(samples) samples[t] = rand(dist[samples[t-1]]) end - Void end """ @@ -314,7 +277,7 @@ The sample path from the `j`-th repetition of the simulation with initial state - `X::Array{Int}` : Array containing the sample paths as columns, of shape (ts_length, k), where k = length(init) if num_reps=nothing, and k = -length(init)*num_reps otherwise. +length(init)* num_reps otherwise. """ function simulate(mc::MarkovChain, @@ -402,7 +365,7 @@ initial states. """ function simulate!(mc::MarkovChain, X::Matrix{Int}) n = size(mc.p, 1) - P_dist = [DiscreteRV(vec(mc.p[i, :])) for i in 1:n] + P_dist = [DiscreteRV(vec(full(mc.p[i, :]))) for i in 1:n] ts_length, k = size(X) diff --git a/src/random_mc.jl b/src/markov/random_mc.jl similarity index 100% rename from src/random_mc.jl rename to src/markov/random_mc.jl diff --git a/test/test_mc_tools.jl b/test/test_mc_tools.jl index d94175a7..a635dac7 100644 --- a/test/test_mc_tools.jl +++ b/test/test_mc_tools.jl @@ -137,35 +137,61 @@ end end @testset "test simulate shape" begin - mc = mc3 - ts_length = 10 - init = [1, 2] - nums_reps = [3, 1] - - @test size(simulate(mc, ts_length)) == (ts_length,) - @test size(simulate(mc, ts_length, init)) == - (ts_length, length(init)) - num_reps = nums_reps[1] - @test size(simulate(mc, ts_length, init, num_reps=num_reps)) == - (ts_length, length(init)*num_reps) - for num_reps in nums_reps + + for mc in (mc3, MarkovChain(sparse(mc3.p))) + ts_length = 10 + init = [1, 2] + nums_reps = [3, 1] + + @test size(simulate(mc, ts_length)) == (ts_length,) + @test size(simulate(mc, ts_length, init)) == + (ts_length, length(init)) + num_reps = nums_reps[1] + @test size(simulate(mc, ts_length, init, num_reps=num_reps)) == + (ts_length, length(init)*num_reps) + for num_reps in nums_reps @test size(simulate(mc, ts_length, num_reps=num_reps)) == - (ts_length, num_reps) + (ts_length, num_reps) + end end end # testset @testset "test simulate init array" begin - mc = mc3 - ts_length = 10 - init = [1, 2] - num_reps = 3 - - X = simulate(mc, ts_length, init, num_reps=num_reps) - @test vec(X[1, :]) == repmat(init, num_reps) + for mc in (mc3, MarkovChain(sparse(mc3.p))) + mc = mc3 + ts_length = 10 + init = [1, 2] + num_reps = 3 + + X = simulate(mc, ts_length, init, num_reps=num_reps) + @test vec(X[1, :]) == repmat(init, num_reps) + end end # testset - @testset "Sparse Matrix " begin + @testset "Sparse Matrix" begin + # use fig1_p from above because we already checked the graph theory + # algorithms for this stochastic matrix. + p = fig1_p + p_s = sparse(p) + @testset "constructors" begin + mc_s = MarkovChain(p_s) + @test isa(mc_s, MarkovChain{eltype(p), typeof(p_s)}) + end + + mc_s = MarkovChain(p_s) + mc = MarkovChain(p) + + @testset "basic correspondence with dense version" begin + @test n_states(mc) == n_states(mc_s) + @test maxabs(mc.p - mc_s.p) == 0.0 + @test recurrent_classes(mc) == recurrent_classes(mc_s) + @test communication_classes(mc) == communication_classes(mc_s) + @test is_irreducible(mc) == is_irreducible(mc_s) + @test is_aperiodic(mc) == is_aperiodic(mc_s) + @test period(mc) == period(mc_s) + @test maxabs(gth_solve(mc_s.p) - gth_solve(mc.p)) < 1e-15 + end end end # testset From 7429bff0b5d1bc23a88b97ae7073532456568b06 Mon Sep 17 00:00:00 2001 From: Spencer Lyon Date: Mon, 14 Mar 2016 09:57:38 -0400 Subject: [PATCH 7/9] make discrete RV fields concrete again --- src/discrete_rv.jl | 12 ++++++------ test/runtests.jl | 4 ---- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/src/discrete_rv.jl b/src/discrete_rv.jl index 09a74c03..a552a9f5 100644 --- a/src/discrete_rv.jl +++ b/src/discrete_rv.jl @@ -28,14 +28,14 @@ vector of probabilities given by q. - `q::Vector{T<:Real}`: A vector of non-negative probabilities that sum to 1 - `Q::Vector{T<:Real}`: The cumulative sum of q """ -type DiscreteRV{T<:Real} - q::AbstractVector{T} - Q::AbstractVector{T} - DiscreteRV(x::AbstractVector{T}) = new(x, cumsum(x)) +type DiscreteRV{TV<:AbstractVector} + q::TV + Q::TV + DiscreteRV(x) = new(x, cumsum(x)) end # outer constructor so people don't have to put {T} themselves -DiscreteRV{T<:Real}(x::AbstractVector{T}) = DiscreteRV{T}(x) +DiscreteRV{TV<:AbstractVector}(x::TV) = DiscreteRV{TV}(x) """ Make a single draw from the discrete distribution @@ -63,4 +63,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] diff --git a/test/runtests.jl b/test/runtests.jl index cd8d1486..5d12c3ff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,10 +34,6 @@ if length(ARGS) > 0 tests = ARGS end -# n = min(8, CPU_CORES, length(tests)) -# n > 1 && addprocs(n) - - srand(42) include("util.jl") From 9ccb9e3d0a40f6f6df27a0eca819fc705caca13e Mon Sep 17 00:00:00 2001 From: Spencer Lyon Date: Wed, 16 Mar 2016 12:21:50 -0400 Subject: [PATCH 8/9] added state_values to MarkovChain also simulate_values and value_simulation --- src/QuantEcon.jl | 6 +- src/discrete_rv.jl | 21 ++-- src/markov/mc_tools.jl | 271 ++++++++++++++++++++--------------------- test/test_mc_tools.jl | 56 ++++++++- 4 files changed, 203 insertions(+), 151 deletions(-) diff --git a/src/QuantEcon.jl b/src/QuantEcon.jl index d22c99cd..026ddc12 100644 --- a/src/QuantEcon.jl +++ b/src/QuantEcon.jl @@ -41,9 +41,11 @@ export # mc_tools MarkovChain, - mc_compute_stationary, mc_sample_path, mc_sample_path!, + mc_compute_stationary, + simulate, simulate!, simulate_values, simulate_values!, + simulation, value_simulation, period, is_irreducible, is_aperiodic, recurrent_classes, - communication_classes, simulate, n_states, + communication_classes, n_states, # gth_solve gth_solve, diff --git a/src/discrete_rv.jl b/src/discrete_rv.jl index a552a9f5..1f30707f 100644 --- a/src/discrete_rv.jl +++ b/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{TV<:AbstractVector} - q::TV - Q::TV - DiscreteRV(x) = 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{TV<:AbstractVector}(x::TV) = DiscreteRV{TV}(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 diff --git a/src/markov/mc_tools.jl b/src/markov/mc_tools.jl index b80a11ba..04b10fa3 100644 --- a/src/markov/mc_tools.jl +++ b/src/markov/mc_tools.jl @@ -27,28 +27,35 @@ transitions. - `p::Matrix` The transition matrix. Must be square, all elements must be positive, and all rows must sum to unity """ -type MarkovChain{T, TM<:AbstractMatrix} +type MarkovChain{T, TM<:AbstractMatrix, TV<:AbstractVector} p::TM # valid stochastic matrix + state_values::TV - function MarkovChain(p::AbstractMatrix) + function MarkovChain(p::AbstractMatrix, state_values) n, m = size(p) - if eltype(p) != T + eltype(p) != T && throw(ArgumentError("Types must be consistent with given matrix")) - elseif n != m - throw(ArgumentError("stochastic matrix must be square")) - elseif any(p.<0) + + n != m && + throw(DimensionMismatch("stochastic matrix must be square")) + + any(p.<0) && throw(ArgumentError("stochastic matrix must have nonnegative elements")) - elseif !check_stochastic_matrix(p) + + !check_stochastic_matrix(p) && throw(ArgumentError("stochastic matrix rows must sum to 1")) - end - return new(p) + length(state_values) != n && + throw(DimensionMismatch("state_values should have $n elements")) + + return new(p, state_values) end end # Provide constructor that infers T from eltype of matrix -MarkovChain(p::AbstractMatrix) = MarkovChain{eltype(p), typeof(p)}(p) +MarkovChain(p::AbstractMatrix, state_values=1:size(p, 1)) = + MarkovChain{eltype(p), typeof(p), typeof(state_values)}(p, state_values) "Number of states in the markov chain `mc`" n_states(mc::MarkovChain) = size(mc.p, 1) @@ -86,7 +93,11 @@ gth_solve{T<:Integer}(a::AbstractMatrix{T}) = # solve x(P-I)=0 by the GTH method function gth_solve{T<:Real}(original::AbstractMatrix{T}) - a = copy(original) + # NOTE: the full here will convert a sparse matrix to dense before the gth + # algorithm begins. Given the nature of the algorithm this will + # almost certainly be more efficient than operating on the sparse + # matrix itself. + a = copy(full(original)) n = size(a, 1) x = zeros(T, n) @@ -176,90 +187,6 @@ function mc_compute_stationary{T}(mc::MarkovChain{T}) return stationary_dists end -""" -Simulate a Markov chain starting from an initial state - -##### Arguments - -- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix -- `init::Int(rand(1:n_states(mc)))` : The index of the initial state. This should -be an integer between 1 and `n_states(mc)` -- `sample_size::Int(1000)`: The number of samples to collect -- `;burn::Int(0)`: The burn in length. Routine drops first `burn` of the -`sample_size` total samples collected - -##### Returns - -- `samples::Vector{Int}`: Vector of simulated states - -""" -function mc_sample_path(mc::MarkovChain, - init::Int=rand(1:n_states(mc)), - sample_size::Int=1000; - burn::Int=0) - samples = Array(Int, sample_size) - samples[1] = init - mc_sample_path!(mc, samples) - samples[burn+1:end] -end - -""" -Simulate a Markov chain starting from an initial distribution - -##### Arguments - -- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix -- `init::Vector` : A vector of length `n_state(mc)` specifying the number -probability of being in seach state in the initial period -- `sample_size::Int(1000)`: The number of samples to collect -- `;burn::Int(0)`: The burn in length. Routine drops first `burn` of the -`sample_size` total samples collected - -##### Returns - -- `samples::Vector{Int}`: Vector of simulated states - -""" -function mc_sample_path(mc::MarkovChain, - init::Vector=vcat(1, zeros(n_states(mc)-1)), - sample_size::Int=1000; burn::Int=0) - # ensure floating point input for Categorical() - init = convert(Vector{Float64}, init) - return mc_sample_path(mc, rand(Categorical(init)), sample_size, burn=burn) -end - -# duplicate python functionality -function mc_sample_path(mc::MarkovChain, - state::Real, - sample_size::Int=1000; burn::Int=0) - init = zeros(n_states(mc)) - init[state] = 1 - return mc_sample_path(mc, init, sample_size, burn=burn) -end - -""" -Fill `samples` with samples from the Markov chain `mc` - -##### Arguments - -- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix -- `samples::Array{Int}` : Pre-allocated vector of integers to be filled with -samples from the markov chain `mc`. The first element will be used as the -initial state and all other elements will be over-written. - -##### Returns - -None modifies `samples` in place -""" -function mc_sample_path!(mc::MarkovChain, samples::Array) - # ensure floating point input for Categorical() - p = convert(Matrix{Float64}, mc.p) - dist = [Categorical(vec(p[i, :])) for i=1:n_states(mc)] - for t=2:length(samples) - samples[t] = rand(dist[samples[t-1]]) - end -end - """ Simulate time series of state transitions of the Markov chain `mc`. @@ -271,24 +198,18 @@ The sample path from the `j`-th repetition of the simulation with initial state - `mc::MarkovChain` : MarkovChain instance. - `ts_length::Int` : Length of each simulation. - `init::Vector{Int}` : Vector containing initial states. -- `;num_reps::Union{Int, Void}(nothing)` : Number of repetitions of simulation. +- `;num_reps::Int(1)` : Number of repetitions of simulation for each element +of `init` ##### Returns -- `X::Array{Int}` : Array containing the sample paths as columns, of shape -(ts_length, k), where k = length(init) if num_reps=nothing, and k = -length(init)* num_reps otherwise. +- `X::Matrix{Int}` : Array containing the sample paths as columns, of shape +(ts_length, k), where k = length(init)* num_reps """ -function simulate(mc::MarkovChain, - ts_length::Int, - init::Vector{Int}; - num_reps::Union{Int, Void}=nothing) - k = length(init) - if num_reps == nothing - num_reps = 1 - end - k *= num_reps +function simulate(mc::MarkovChain, ts_length::Int, init::Vector{Int}; + num_reps::Int=1) + k = length(init)*num_reps X = Array(Int, ts_length, k) X[1, :] = repmat(init, num_reps) @@ -304,26 +225,16 @@ Simulate time series of state transitions of the Markov chain `mc`. - `mc::MarkovChain` : MarkovChain instance. - `ts_length::Int` : Length of each simulation. - `init::Int` : Initial state. -- `;num_reps::Union{Int, Void}(nothing)` : Number of repetitions of simulation. +- `;num_reps::Int(1)` : Number of repetitions of simulation ##### Returns -- `X::Array{Int}` : Array containing the sample path(s) as column(s), of shape -(ts_length,) if num_reps=nothing; of shape (ts_length, num_reps) otherwise. +- `X::Matrix{Int}` : Array containing the sample paths as columns, of shape +(ts_length, k), where k = num_reps """ -function simulate(mc::MarkovChain, - ts_length::Int, - init::Int; - num_reps::Union{Int, Void}=nothing) - init_vec = [init] - X = simulate(mc, ts_length, init_vec, num_reps=num_reps) - if num_reps == nothing - X = vec(X) - end - - return X -end +simulate(mc::MarkovChain, ts_length::Int, init::Int; num_reps::Int=1) = + simulate(mc, ts_length, [init], num_reps=num_reps) """ Simulate time series of state transitions of the Markov chain `mc`. @@ -336,19 +247,12 @@ Simulate time series of state transitions of the Markov chain `mc`. ##### Returns -- `X::Array{Int}` : Array containing the sample path(s) as column(s), of shape -(ts_length,) if num_reps=nothing; of shape (ts_length, num_reps) otherwise, -where the initial state is randomly drawn for each simulation. +- `X::Matrix{Int}` : Array containing the sample paths as columns, of shape +(ts_length, k), where k = num_reps """ -function simulate(mc::MarkovChain, - ts_length::Int; - num_reps::Union{Int, Void}=nothing) - if num_reps == nothing - init = rand(1:size(mc.p, 1)) - else - init = rand(1:size(mc.p, 1), num_reps) - end +function simulate(mc::MarkovChain, ts_length::Int; num_reps::Int=1) + init = rand(1:size(mc.p, 1), num_reps) return simulate(mc, ts_length, init) end @@ -359,13 +263,17 @@ Fill `X` with sample paths of the Markov chain `mc` as columns. - `mc::MarkovChain` : MarkovChain instance. - `X::Matrix{Int}` : Preallocated matrix of integers to be filled with sample -paths of the markov chain `mc`. The elements in `X[0, :]` will be used as the +paths of the markov chain `mc`. The elements in `X[1, :]` will be used as the initial states. """ function simulate!(mc::MarkovChain, X::Matrix{Int}) n = size(mc.p, 1) - P_dist = [DiscreteRV(vec(full(mc.p[i, :]))) for i in 1:n] + + # NOTE: ensure dense array and transpose before slicing the array. Then + # when passing to DiscreteRV use `sub` to avoid allocating again + p = full(mc.p)' + P_dist = [DiscreteRV(sub(p, :, i)) for i in 1:n] ts_length, k = size(X) @@ -374,4 +282,95 @@ function simulate!(mc::MarkovChain, X::Matrix{Int}) X[t+1, i] = draw(P_dist[X[t]]) end end + X +end + +""" +Simulate time series of state transitions of the Markov chain `mc`. + +##### Arguments + +- `mc::MarkovChain` : MarkovChain instance. +- `ts_length::Int` : Length of each simulation. +- `init_state::Int(rand(1:n_states(mc)))` : Initial state. + +##### Returns + +- `x::Vector`: A vector of transition indices for a single simulation +""" +function simulation(mc::MarkovChain, ts_length::Int, + init_state::Int=rand(1:n_states(mc))) + simulate(mc, ts_length, init_state; num_reps=1)[:, 1] +end + +# ----------------------- # +# simulate_values methods # +# ----------------------- # + +function simulate_values(mc::MarkovChain, + ts_length::Int, + init_state::Vector{Int}; + num_reps::Int=1) + k = length(init_state)*num_reps + X = Array(eltype(mc.state_values), ts_length, k) + init_state = repmat(init_state, num_reps) + X[1, :] = mc.state_values[init_state] + + simulate_values!(mc, init_state, X) + return X +end + +function simulate_values(mc::MarkovChain, ts_length::Int, init_state::Int; + num_reps::Int=1) + simulate_values(mc, ts_length, [init_state], num_reps=num_reps) +end + +function simulate_values(mc::MarkovChain, ts_length::Int; num_reps::Int=1) + init_state = rand(1:size(mc.p, 1), num_reps) + simulate_values(mc, ts_length, init_state) +end + +function simulate_values!(mc::MarkovChain, init_state::Vector{Int}, X::Matrix) + n = size(mc.p, 1) + + # NOTE: ensure dense array and transpose before slicing the array. Then + # when passing to DiscreteRV use `sub` to avoid allocating again + p = full(mc.p)' + P_dist = [DiscreteRV(sub(p, :, i)) for i in 1:n] + + ts_length, k = size(X) + + for i in 1:k + i_state = init_state[i] + for t in 1:ts_length-1 + i_state = draw(P_dist[i_state]) + X[t+1, i] = mc.state_values[i_state] + end + end + X +end + +@doc """ Like `simulate(::MarkovChain, args...; kwargs...)`, but instead of +returning integers specifying the state indices, this routine returns the +values of the `mc.state_values` at each of those indices. See docstring +for `simulate` for more information +""" simulate_values, simulate_values! + + +""" +Simulate time series of state transitions of the Markov chain `mc`. + +##### Arguments + +- `mc::MarkovChain` : MarkovChain instance. +- `ts_length::Int` : Length of each simulation. +- `init_state::Int(rand(1:n_states(mc)))` : Initial state. + +##### Returns + +- `x::Vector`: A vector of state values along a simulated path +""" +function value_simulation(mc::MarkovChain, ts_length::Int, + init_state::Int=rand(1:n_states(mc))) + simulate_values(mc, ts_length, init_state; num_reps=1)[:, 1] end diff --git a/test/test_mc_tools.jl b/test/test_mc_tools.jl index a635dac7..b0e32c86 100644 --- a/test/test_mc_tools.jl +++ b/test/test_mc_tools.jl @@ -113,9 +113,17 @@ end end @testset "test MarkovChain throws errors" begin - @test_throws ArgumentError MarkovChain(rand(4, 5)) # not square - @test_throws ArgumentError MarkovChain([0.0 0.5; 0.2 0.8]) # first row doesn't sum to 1 - @test_throws ArgumentError MarkovChain([-1 1; 0.2 0.8]) # negative element, but sums to 1 + # not square + @test_throws DimensionMismatch MarkovChain(rand(4, 5)) + + # state_values to long + @test_throws DimensionMismatch MarkovChain([0.5 0.5; 0.2 0.8], rand(3)) + + # first row doesn't sum to 1 + @test_throws ArgumentError MarkovChain([0.0 0.5; 0.2 0.8]) + + # negative element, but sums to 1 + @test_throws ArgumentError MarkovChain([-1 1; 0.2 0.8]) end @testset "test graph theoretic algorithms" begin @@ -143,7 +151,7 @@ end init = [1, 2] nums_reps = [3, 1] - @test size(simulate(mc, ts_length)) == (ts_length,) + @test size(simulate(mc, ts_length)) == (ts_length,1) @test size(simulate(mc, ts_length, init)) == (ts_length, length(init)) num_reps = nums_reps[1] @@ -158,7 +166,6 @@ end @testset "test simulate init array" begin for mc in (mc3, MarkovChain(sparse(mc3.p))) - mc = mc3 ts_length = 10 init = [1, 2] num_reps = 3 @@ -194,4 +201,43 @@ end end end + @testset "simulate_values" begin + for _mc in (mc3, MarkovChain(sparse(mc3.p))) + for T in [Float16, Float32, Float64, Int8, Int16, Int32, Int64] + mc = MarkovChain(_mc.p, rand(T, size(_mc.p, 1))) + ts_length = 10 + init = [1, 2] + nums_reps = [3, 1] + + @test size(@inferred(simulate_values(mc, ts_length))) == (ts_length,1) + @test size(@inferred(simulate_values(mc, ts_length, init))) == + (ts_length, length(init)) + num_reps = nums_reps[1] + @test size(simulate(mc, ts_length, init, num_reps=num_reps)) == + (ts_length, length(init)*num_reps) + for num_reps in nums_reps + @test size(simulate(mc, ts_length, num_reps=num_reps)) == + (ts_length, num_reps) + end + end # state_values eltypes + end + + @testset "test simulate_values init array" begin + for mc in (mc3, MarkovChain(sparse(mc3.p))) + ts_length = 10 + init = [1, 2] + num_reps = 3 + + X = simulate(mc, ts_length, init, num_reps=num_reps) + @test vec(X[1, :]) == repmat(init, num_reps) + end + end # testset + end + + @testset "simulation and value_simulation" begin + @test size(@inferred(simulation(mc3, 10, 1))) == (10,) + @test size(@inferred(value_simulation(mc3, 10, 1))) == (10,) + end + + end # testset From 44a7314830556281e258b25d48a9e4814cad1d7f Mon Sep 17 00:00:00 2001 From: Spencer Lyon Date: Wed, 16 Mar 2016 13:56:20 -0400 Subject: [PATCH 9/9] tauchen and rouwenhorst now return MarkovChain instances --- src/markov/markov_approx.jl | 23 ++++++++++++++--------- test/test_markov_approx.jl | 21 +++++---------------- 2 files changed, 19 insertions(+), 25 deletions(-) diff --git a/src/markov/markov_approx.jl b/src/markov/markov_approx.jl index 65ebcb90..cdf0f469 100644 --- a/src/markov/markov_approx.jl +++ b/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] + diff --git a/test/test_markov_approx.jl b/test/test_markov_approx.jl index 01cc0a00..43bb133f 100644 --- a/test/test_markov_approx.jl +++ b/test/test_markov_approx.jl @@ -5,26 +5,15 @@ ρ, σ_u = rand(2) μ = 0.0 - + for n=3:25 - x, P = rouwenhorst(n, ρ, σ_u, μ) + mc = rouwenhorst(n, ρ, σ_u, μ) - @test size(P, 1) == size(P, 2) - @test ndims(x) == 1 - @test ndims(P) == 2 - @test isapprox(sum(P, 2), ones(n, 1); rough_kwargs...) - @test all(P .>= 0.0) == true - @test isapprox(sum(x) , 0.0; rough_kwargs...) + @test isapprox(sum(mc.state_values) , 0.0; rough_kwargs...) for m=1:4 - x, P = tauchen(n, ρ, σ_u, μ, m) - - @test size(P, 1) == size(P, 2) - @test ndims(x) == 1 - @test ndims(P) == 2 - @test isapprox(sum(P, 2), ones(n, 1); rough_kwargs...) - @test all(P .>= 0.0) == true - @test isapprox(sum(x), 0.0; rough_kwargs...) + mc = tauchen(n, ρ, σ_u, μ, m) + @test isapprox(sum(mc.state_values), 0.0; rough_kwargs...) end end end # @testset