Violante HW 3 Code

readme.md

This gist contains the code used to define the utilities needed to complete homework 3 for Gianluca Violante's first year macro course and NYU.

discretize.jl

"""
Various routines to discretize AR(1) processes

@author : Spencer Lyon <spencer.lyon@nyu.edu>
@date : 2014-04-10 23:55:05

"""
norm_cdf{T <: Real}(x::Union(T, Array{T})) = 0.5 * erfc(-x/sqrt(2))


function tauchen(N::Int64, ρ::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::Int, ρ::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

driver.jl

# Run this file by calling the command `julia -p 7 -L endog_grid driver`

standard_mod = IFP()
models = {standard_mod,  # Part a (and b/c) -- standard parameters
          IFP(gamma=1.0),  # part b
          IFP(gamma=5.0),  # part b
          IFP(sig_eps=0.01),  # part c
          IFP(sig_eps=0.12),  # part c
          IFP(a_min=-minimum(standard_mod.y)/standard_mod.r),  # part d
          }

# solve each model in parallel
models = pmap(endog_grid!, models)

# simulate each model in parallel
@everywhere sim_pan(m::IFP) = simulation(m, N=250, T=500, burn=99);
sims = pmap(sim_pan, models)

endog_grid.jl

"""
Solving the consumption savings (income fluctuation) problem using the
endogenous grid method proposed by Carroll (2005)

@author : Spencer Lyon <spencer.lyon@nyu.edu>
@date : 2014-04-14 17:38:35

"""
using Grid
require("patch.jl")
require("ifp.jl")

# Define utility function and marginal utility and its inverse
# NOTE: These will all be inlined
u{T <: Real}(c::Union(T, Array{T}), γ::Float64) = c.^(1.0 - γ) ./ (1.0 - γ)
up{T <: Real}(c::Union(T, Array{T}), γ::Float64) = c.^(- γ)
up_inv{T <: Real}(b::Union(T, Array{T}), γ::Float64) = b.^(- 1.0 / γ)


#=
    Set up the initial guess for the policy function of c.

    It uses the analytical solution to the same problem with quadratic
    utility and iid shocks.
=#
function initialize_c(a::Vector{Float64}, y::Vector{Float64}, r::Float64)
    n_a, n_y = map(length, {a, y})
    c = Array(Float64, n_a, n_y)
    for i=1:n_a, j=1:n_y
        c[i, j] = r*a[i] + y[j]
    end
    return c
end


#=
    Update the guess for c using linear interpolation. This routine is
    called at the end of each loop in the main endogenous grid
    algorithm.
=#
# Updates c_til_i inplace
function update_c!(a_star::Matrix{Float64}, c_til::Matrix{Float64},
                   c_til_i::Matrix{Float64}, a_star_1::Vector{Float64},
                   a::Vector{Float64}, y::Vector{Float64}, R::Float64)

    num_bind::Int64 = 0
    n_a, n_y = size(c_til)
    for j=1:n_y  # Do interpolation (step 7)
        ig = InterpIrregular(a_star[:, j], c_til[:, j], BCnearest,
                             InterpLinear)
        c_til_i[:, j] = ig[a]

        # update boundary condition on lower end
        for i=1:n_a
            if a[i] < a_star_1[j]
                num_bind += 1
                c_til_i[i, j] = R * a[i] + y[j] - a[1]
            end
        end
    end
    return num_bind
end


#=
    Apply the endogenous grid method to solve for the policy functions
    c(a, y) and a'(a, y) for a given specification of the coefficient of
    absolute risk aversion. The return value is two (n_a, n_y) matrices
    representing the optimal policy rules for each value on the grids
    for a and y.

    tol tolerance level for convergence
=#
function endog_grid!(model::IFP; tol::Float64=1e-12, verbose=false)
    # Unpack parameters from model object
    n_a::Integer = model.n_a
    n_y::Integer = model.n_y
    a::Vector{Float64} = model.a

    # step 1 (discretize process for y) was done when model obj was created
    y::Vector{Float64} = model.y
    P::Matrix{Float64} = model.P
    beta::Float64 = model.beta
    gamma::Float64 = model.gamma
    r::Float64 = model.r
    R::Float64 = model.R

    # Allocate memory for arrays used in iteration
    B = Array(Float64, n_a, n_y)
    c_til = Array(Float64, n_a, n_y)
    a_star = Array(Float64, n_a, n_y)
    a_star_1 = Array(Float64, n_y)
    c_til_i = Array(Float64, n_a, n_y)
    c = initialize_c(a, y, r)  # step 2

    # set up iteration parameters
    err = 1.0
    it = 0

    # Main algorithm
    while err > tol
        B[:] = beta * R .* (up(c, gamma) * P')  # step 3
        c_til[:] = up_inv(B, gamma)::Matrix{Float64}  # step 4
        a_star[:] = (c_til .+ a .- y') ./ R  # step 5
        a_star_1[:] = a_star[1, :]  # step 6
        bound = update_c!(a_star, c_til, c_til_i, a_star_1, a, y, R)  # step 7
        err = maximum(abs(c - c_til_i))  # step 8

        # Update guess for c
        c[:] = c_til_i

        # print status, break if running too long
        it += 1
        if verbose && it % 5 == 0
            println("it=$it\tbound=$bound\terr=$err")
        end
        if it > 300 break end
    end

    println("Converged. Total iterations: $it")

    # compute a'(a, y) from budget constraint using c
    ap = (R .* a) .+  y' - c

    # Attach policy rules to model object
    model.Sol = IFPSol(c, ap)

    return model # don't return anything. Just update model inplace
end


#=
    Given a set of policy functions, the Markov process for y,
    and a desired simulation length, simulate the economy for a
    specified number of periods. The return value is three Vectors: one
    for the path of consumption, one for the path of asset holdings,
    and a third for the path of the endowment process.

    a0_i is the index of the starting level of asset holdings.

    burn is the burn in period. It is set to remove 10% of simulation
    length by default.

    seed is the seed for the random number generator. The default value
    is a random number (what?!! bootstrapping the RNG?!)
=#
function simulation(model::IFP; T::Integer=20000, a0::Number=model.a_min,
                    N::Integer=1, burn::Integer=2000,
                    seed::Number=int(500*rand()), use_df=false)
    if ~isdefined(model, :Sol)
        msg = "Sol field not found on model.\n"
        msg *= "Must call endog_grid! before simulation"
        msg *= "I'll try to do it for you"
        println(msg)
        endog_grid!(model, verbose=false)
    end

    # Pull out needed objects from model
    c::Matrix{Float64} = model.Sol.c
    a_prime::Matrix{Float64} = model.Sol.ap
    m::Markov = model.m
    n_a::Integer = model.n_a
    n_y::Integer = model.n_y
    R::Float64 = model.R

    # Set up data storage
    c_data = Array(Float64, T, N)
    y_data = Array(Float64, T, N)
    a_data = Array(Float64, T+1, N)

    # seed RNG to get consistent simulation of Markov process
    srand(seed)

    # Construct interpolator objects to use in computing c inside loop
    igs = [InterpIrregular(a_prime[:, j], c[:, j], BCnearest, InterpLinear)
           for j=1:n_y]

    i_t = Array(Int8, T)
    for i=1:N
        # simulate markov process
        y_t, i_t = simulate(m, T, ret_ind=true)
        y_data[:, i] = y_t

        # use initial condition to start income process
        a_data[1, i] = a0

        # simulate forward all T time periods
        for t=1:T
            c_data[t, i] = igs[i_t[t]][a_data[t, i]]
            a_data[t+1, i] = R * a_data[t, i] + y_data[t, i] - c_data[t, i]
        end
    end

    # if use_df  # Convert to DataFrame at end
    #     c_data = convert(DataFrame, c_data)
    #     names!(c_data, [symbol("csim_$i") for i=1:50])
    #     y_data = convert(DataFrame, y_data)
    #     names!(y_data, [symbol("ysim_$i") for i=1:50])
    #     a_data = convert(DataFrame, a_data)
    #     names!(a_data, [symbol("asim_$i") for i=1:50])
    # end

    # burn = min(burn, int(0.1*T))  # Don't burn more than we have
    return c_data[burn+1:end, :], a_data[burn+1:end, :], y_data[burn+1:end, :]

end

hw3.jl

"""
Computational homework assignment

@author : Spencer Lyon <spencer.lyon@nyu.edu>
@date : 2014-04-10 23:55:05

"""
require("discretize.jl")
require("markov.jl")

# Useful functions
logn_E(μ, σ) = μ + σ^2/2
logn_var(μ, σ) = (exp(σ^2) - 1) * exp(2μ + σ^2)

function sim_moments(m::Markov; T=50000, burn=2000)
    """
    Given an object of type Markov, simulate a markov process, compute
    mean and variance of simulation, and compare to theoretical moments.
    """
    # Simulate and compute moments (throw out burn in period)

    sim = simulate(m, T)

    mu = mean(sim[burn:end])
    sig2 = var(sim[burn:end])

    # print nice string representation of moments
    E_str = "\t" * rpad(@sprintf("E[y]   = %.6f", mu), 17)
    E_str *= "\t" * @sprintf("Expected = %.6f", μ_y)
    E_str *= "\t" * @sprintf("Abs. Difference = %.6f", abs(mu - μ_y))
    V_str = "\n\t" * rpad(@sprintf("Var[y] = %.6f", sig2), 17)
    V_str *= "\t" * @sprintf("Expected = %.6f", σ2_y)
    V_str *= "\t" * @sprintf("Abs. Difference = %.6f", abs(sig2 - σ2_y))

    println(E_str * V_str)

    return mu, sig2, sim
end

type Result
    w
    P
    mu
    sig2
end

# Define parameters
ρ = 0.90
σ_ε = 0.06
wbar = - σ_ε^2 / (2 * (1 + ρ))

# theoretical moments
μ_w = wbar / (1 - ρ)
σ_w = sqrt(σ_ε ^2 / (1 - ρ^2))
μ_y = 1.0
σ2_y = logn_var(μ_w, σ_w)

## -------------------------- ##
#- Parts a, b: Tauchen Method -#
## -------------------------- ##

t_res = Dict{Int, Result}()
for N in [5, 10]
    # Use tauchen method to estimate AR(1) with markov process
    w, P = tauchen(N, ρ, σ_ε, wbar, 3)
    println("Tauchen method with $N state Markov process:")
    mu, sig2, sim = sim_moments(Markov(P, exp(w)))

    # Add Result object to dict
    t_res[N] = Result(w, P, mu, sig2)
end

## ---------------------------------- ##
#- Part c: 5 state Rouwenhorst Method -#
## ---------------------------------- ##

# Use rouwenhorst to approximate process
w_c, P_c = rouwenhorst(5, ρ, σ_ε, wbar)
println("Rouwenhorst method with 5 state Markov process:")
mu_c, sig2_c, sim_c = sim_moments(Markov(P_c, exp(w_c)))

## ----------------------------------- ##
#- Part d: Rouwenhorst, new parameters -#
## ----------------------------------- ##

# Define new parameters
ρ_d = 0.98
σ_ε_d = σ_ε * sqrt(1 - ρ_d^2) / sqrt(1 - ρ^2)
wbar_d = - σ_ε_d^2 / (2 * (1 + ρ_d))

w_d, P_d = rouwenhorst(5, ρ_d, σ_ε_d, wbar_d)
println("Rouwenhorst method with 5 state Markov process (new params):")
mu_d, sig2_d, sim_d = sim_moments(Markov(P_d, exp(w_d)))

ifp.jl

"""
Type to define the income fluctuation problem (consumption-savings
problem)

@author : Spencer Lyon <spencer.lyon@nyu.edu>
@date : 2014-04-19 20:10:42

"""
require("markov.jl")
require("discretize.jl")

import Base.show

abstract Model
abstract Solution

type IFPSol <: Solution
    c::Matrix{Float64}
    ap::Matrix{Float64}

    function IFPSol(c::Matrix{Float64}, ap::Matrix{Float64})
        new(c, ap)
    end
end


type IFP{T <: FloatingPoint} <: Model
    # Model parameters
    rho::T
    beta::T
    r::T
    R::T
    gamma::T

    # Grid parameters for a
    n_a::Integer
    a_min::T
    a_max::T
    a::Vector{T}

    # Parameters for y
    n_y::Integer
    sig_eps::T
    y::Vector{T}
    P::Matrix{T}
    m::Markov

    # Hold the solution
    Sol::IFPSol

    # Incomplete initialization. Everything but Sol field
    function IFP(rho::T, beta::T, r::T, R::T, gamma::T,
                 n_a::Integer, a_min::T, a_max::T, a::Vector{T},
                 n_y::Integer, sig_eps::T, y::Vector{T}, P::Matrix{T},
                 m::Markov)

        new(rho, beta, r, R, gamma, n_a, a_min, a_max, a,
            n_y, sig_eps, y, P, m)
    end

    # Default constructor
    function IFP(rho::T, beta::T, r::T, R::T, gamma::T,
                 n_a::Integer, a_min::T, a_max::T, a::Vector{T},
                 n_y::Integer, sig_eps::T, y::Vector{T}, P::Matrix{T},
                 m::Markov, Sol::IFPSol)

        new(rho, beta, r, R, gamma, n_a, a_min, a_max, a, n_y,
            sig_eps, y, P, m, Sol)
    end

end


function IFP{T <: FloatingPoint}(;rho::T=0.90, beta::T=0.95, r::T=0.02,
                                 gamma::T=2.0,
                                 n_a::Integer=300, a_min::T=0.0, a_max::T=100.,
                                 n_y::Integer=5, sig_eps::T=0.06)

    # Compute a
    a::Vector{T} = linspace(a_min, a_max, n_a)

    # Compute R
    R = 1.0 + r

    # compute P, y, m
    wbar = - sig_eps^2 / (2 * (1 + rho))
    w::Vector{T}, P::Matrix{T} = rouwenhorst(n_y, rho, sig_eps, wbar)
    y::Vector{T} = exp(w)
    m::Markov = Markov(P, y)

    return IFP{T}(rho, beta, r, R, gamma, n_a, a_min, a_max, a, n_y, sig_eps,
                  y, P, m)
end


function show(io::IO, ifp::IFP)
    msg = "Income Fluctuation Problem\n"
    msg *= "-" ^ (length(msg) - 1)

    param_str = ""
    param_syms = names(ifp)
    n = maximum(map(length, map(string, param_syms)))
    for p in param_syms
        if is(p, :Sol)
            if isdefined(ifp, :Sol)
                param_str *= "\t$(rpad("Sol", n)): Computed\n"
            else
                param_str *= "\t$(rpad("Sol", n)): Not yet computed\n"
            end
            continue
        end

        f = getfield(ifp, p)
        if isa(f, Array)
            param_str *= "\t$(rpad(string(p), n)): $(summary(f))\n"
        else
            param_str *= "\t$(rpad(string(p), n)): $f\n"
        end

    end

    msg *= "\n" * param_str

    print(io, msg)
end

markov.jl

import Base.show

type Markov
    P::Matrix{Float64}
    y::Vector{Float64}   # vector of state values
    pi_0::Union(Vector{Float64}, Vector{Int64})

    function Markov(P, y, pi_0)
        if size(P, 1) != size(P, 2)
            error("P must be a square matrix")
        end

        if size(pi_0, 1) != size(P, 1)
            error("pi_0 must have same number of elements as each row of P")
        end

        if size(y, 1) != size(P, 1)
            error("y must have same number of elements as each row of P")
        end

        new(P, y, pi_0)
    end
end


function Markov(P::Matrix{Float64}, y::Vector{Float64})
    pi_0 = zeros(size(P, 1))
    pi_0[1] = 1.0
    Markov(P, y, pi_0)
end


function Markov(P::Matrix{Float64})
    y = [1.:size(P, 1)]
    pi_0 = zeros(size(P, 1))
    pi_0[1] = 1.0
    Markov(P, y, pi_0)
end


function show(io::IO, m::Markov)
    msg = "$(length(m.y)) state Markov Chain"
    print(io, msg)
end


function simulate(mc::Markov, N::Int; ret_ind=false)
    # take cumsum along each row of transition matrix
    cs_P = cumsum(mc.P, 2)

    # Draw N random numbers on U[0, 1) and allocate memory for draw
    α = rand(N)
    draw = Array(Int16, N)

    # Compute starting value using initial distribution pi_0
    draw[1] = findfirst(α[1] .< cumsum(mc.pi_0))

    # Fill in chain using transition matrix
    for j=2:N
        draw[j] = findfirst(α[j] .< cs_P[draw[j-1], :])
    end

    if ret_ind
        return mc.y[draw], draw
    else
        return mc.y[draw]
    end
end

patch.jl

import Base.getindex

using Grid

## Vectorized evaluation at multiple points
function getindex{T,R<:Real}(G::InterpIrregular{T,1}, x::AbstractVector{R})
    n = length(x)
    v = Array(T, n)
    for i = 1:n
        v[i] = getindex(G, x[i])
    end
    v
end