ENH: moved tauchen routine to mc_tools. updated tests

src/QuantEcon.jl

@@ -41,14 +41,8 @@ export
 
 # mc_tools
     MarkovChain,
-    mc_compute_stationary, mc_sample_path, mc_sample_path!,
-
-# gth_solve
-    gth_solve,
-
-# markov_approx
-    tauchen,
-    rouwenhorst,
+    mc_compute_stationary, mc_sample_path, mc_sample_path!, gth_solve,
+    tauchen, rouwenhorst,
 
 # lae
     LAE,
@@ -116,7 +110,6 @@ include("lae.jl")
 include("lqcontrol.jl")
 include("lqnash.jl")
 include("lss.jl")
-include("markov_approx.jl")
 include("matrix_eqn.jl")
 include("mc_tools.jl")
 include("robustlq.jl")

src/mc_tools.jl

@@ -191,4 +191,141 @@ function mc_sample_path!(mc::MarkovChain, samples::Vector)
     length(samples) < 2 &&
         throw(ArgumentError("samples vector must have length greater than 2"))
     samples = mc_sample_path(mc,samples[1],length(samples)-1)
-end
+end
+
+## ----------------------------- ##
+#- AR(1) discretization routines -#
+## ----------------------------- ##
+
+norm_cdf{T <: Real}(x::T) = 0.5 * erfc(-x/sqrt(2))
+norm_cdf{T <: Real}(x::Array{T}) = 0.5 .* erfc(-x./sqrt(2))
+
+
+function tauchen(N::Integer, ρ::Real, σ::Real, μ::Real=0.0, n_std::Int64=3)
+    """
+    Use Tauchen's (1986) method to produce finite state Markov
+    approximation of the AR(1) processes
+
+    y_t = μ + ρ y_{t-1} + ε_t,
+
+    where ε_t ~ N (0, σ^2)
+
+    Parameters
+    ----------
+    N : int
+        Number of points in markov process
+
+    ρ : float
+        Persistence parameter in AR(1) process
+
+    σ : float
+        Standard deviation of random component of AR(1) process
+
+    μ : float, optional(default=0.0)
+        Mean of AR(1) process
+
+    n_std : int, optional(default=3)
+        The number of standard deviations to each side the processes
+        should span
+
+    Returns
+    -------
+    y : array(dtype=float, ndim=1)
+        1d-Array of nodes in the state space
+
+    Π : array(dtype=float, ndim=2)
+        Matrix transition probabilities for Markov Process
+
+    """
+    # Get discretized space
+    a_bar = n_std * sqrt(σ^2 / (1 - ρ^2))
+    y = linspace(-a_bar, a_bar, N)
+    d = y[2] - y[1]
+
+    # Get transition probabilities
+    Π = zeros(N, N)
+    for row = 1:N
+        # Do end points first
+        Π[row, 1] = norm_cdf((y[1] - ρ*y[row] + d/2) / σ)
+        Π[row, N] = 1 - norm_cdf((y[N] - ρ*y[row] - d/2) / σ)
+
+        # fill in the middle columns
+        for col = 2:N-1
+            Π[row, col] = (norm_cdf((y[col] - ρ*y[row] + d/2) / σ) -
+                           norm_cdf((y[col] - ρ*y[row] - d/2) / σ))
+        end
+    end
+
+    # NOTE: I need to shift this vector after finding probabilities
+    #       because when finding the probabilities I use a function norm_cdf
+    #       that assumes its input argument is distributed N(0, 1). After
+    #       adding the mean E[y] is no longer 0, so I would be passing
+    #       elements with the wrong distribution.
+    #
+    #       It is ok to do after the fact because adding this constant to each
+    #       term effectively shifts the entire distribution. Because the
+    #       normal distribution is symmetric and we just care about relative
+    #       distances between points, the probabilities will be the same.
+    #
+    #       I could have shifted it before, but then I would need to evaluate
+    #       the cdf with a function that allows the distribution of input
+    #       arguments to be [μ/(1 - ρ), 1] instead of [0, 1]
+    y .+= μ / (1 - ρ) # center process around its mean (wbar / (1 - rho))
+
+    return y, Π
+end
+
+
+function rouwenhorst(N::Integer, ρ::Real, σ::Real, μ::Real=0.0)
+    """
+    Use Rouwenhorst's method to produce finite state Markov
+    approximation of the AR(1) processes
+
+    y_t = μ + ρ y_{t-1} + ε_t,
+
+    where ε_t ~ N (0, σ^2)
+
+    Parameters
+    ----------
+    N : int
+        Number of points in markov process
+
+    ρ : float
+        Persistence parameter in AR(1) process
+
+    σ : float
+        Standard deviation of random component of AR(1) process
+
+    μ : float, optional(default=0.0)
+        Mean of AR(1) process
+
+    Returns
+    -------
+    y : array(dtype=float, ndim=1)
+        1d-Array of nodes in the state space
+
+    Θ : array(dtype=float, ndim=2)
+        Matrix transition probabilities for Markov Process
+
+    """
+    σ_y = σ / sqrt(1-ρ^2)
+    p  = (1+ρ)/2
+    Θ = [p 1-p; 1-p p]
+
+    for n = 3:N
+        z_vec = zeros(n-1,1)
+        z_vec_long = zeros(1, n)
+        Θ = p.*[Θ z_vec; z_vec_long] +
+            (1-p).*[z_vec Θ; z_vec_long] +
+            (1-p).*[z_vec_long; Θ z_vec] +
+            p.*[z_vec_long; z_vec Θ]
+        Θ[2:end-1,:] ./=  2.0
+    end
+
+    ψ = sqrt(N-1) * σ_y
+    w = linspace(-ψ, ψ, N)
+    w .+= μ / (1 - ρ)  # center process around its mean (wbar / (1 - rho))
+    return w, Θ
+end
+
+

test/runtests.jl

@@ -12,7 +12,6 @@ include("test_lae.jl")
 include("test_lqcontrol.jl")
 include("test_lqnash.jl")
 include("test_lss.jl")
-include("test_markov_approx.jl")
 include("test_matrix_eqn.jl")
 include("test_mc_tools.jl")
 include("test_models.jl")

test/test_mc_tools.jl

@@ -102,6 +102,24 @@ facts("Testing mc_tools.jl") do
         @fact_throws MarkovChain([0.0 0.5; 0.2 0.8])  # first row doesn't sum to 1
         @fact_throws MarkovChain([-1 1; 0.2 0.8])  # negative element, but sums to 1
     end
-end  # facts
 
-end  # module
+    context("testing tauchen method") do
+    for n=3:25, m=1:4
+        rough_kwargs = {:atol => 1e-8, :rtol => 1e-8}
+
+        # set up
+        ρ, σ_u = rand(2)
+        μ = 0.0
+        x, P = tauchen(n, ρ, σ_u, μ, m)
+
+        @fact size(P, 1) => size(P, 2)
+        @fact ndims(x) => 1
+        @fact ndims(P) => 2
+        @fact sum(P, 2) => roughly(ones(n, 1); rough_kwargs...)
+        @fact all(P .>= 0.0) => true
+        @fact sum(x) => roughly(0.0; rough_kwargs...)
+        end
+    end  # context
+
+end  # facts
+end  # module