Merge branch 'master' of https://github.com/QuantEcon/QuantEcon.jl

REQUIRE

@@ -5,5 +5,5 @@ LightGraphs
 Primes
 Compat 0.18.0
 StatsBase
-Optim 0.7.3
+Optim 0.9.0
 NLopt 0.3.4

src/QuantEcon.jl

@@ -53,6 +53,9 @@ export
     communication_classes, n_states,
     discrete_var, Even,
 
+# modeltools
+    AbstractUtility, LogUtility, CRRAUtility, CFEUtility, EllipticalUtility, derivative,
+
 # gth_solve
     gth_solve,
 
@@ -156,5 +159,6 @@ include("quad.jl")
 include("quadsums.jl")
 include("zeros.jl")
 include("optimization.jl")
+include("modeltools/utility.jl")
 
 end # module

src/markov/markov_approx.jl

@@ -616,8 +616,8 @@ function discrete_approximation(D::AbstractVector, T::Function, TBar::AbstractVe
     # Compute maximum entropy discrete distribution
     options = Optim.Options(f_tol=1e-16, x_tol=1e-16)
     obj(lambda) = entropy_obj(lambda, Tx, TBar, q)
-    grad!(lambda, grad) = entropy_grad!(grad, lambda, Tx, TBar, q)
-    hess!(lambda, hess) = entropy_hess!(hess, lambda, Tx, TBar, q)
+    grad!(grad, lambda) = entropy_grad!(grad, lambda, Tx, TBar, q)
+    hess!(hess, lambda) = entropy_hess!(hess, lambda, Tx, TBar, q)
     res = Optim.optimize(obj, grad!, hess!, lambda0, Optim.Newton(), options)
     # Sometimes the algorithm fails to converge if the initial guess is too far
     # away from the truth. If this occurs, the program tries an initial guess
@@ -635,7 +635,7 @@ function discrete_approximation(D::AbstractVector, T::Function, TBar::AbstractVe
     Tdiff = Tx .- TBar
     p = (q'.*exp.(lambda_bar'*Tdiff))/minimum_value
     grad = similar(lambda0)
-    grad!(Optim.minimizer(res), grad)
+    grad!(grad, Optim.minimizer(res))
     moment_error = grad/minimum_value
     return p, lambda_bar, moment_error
 end

src/modeltools/utility.jl

@@ -0,0 +1,110 @@
+abstract type AbstractUtility end
+
+#
+# Separable utility
+#
+
+# Consumption utility
+
+"""
+Type used to evaluate log utility. Log utility takes the form
+
+u(c) = \log(c)
+
+Additionally, this code assumes that if c < 1e-10 then
+
+u(c) = log(1e-10) + 1e10*(c - 1e-10)
+
+"""
+struct LogUtility <: AbstractUtility
+    ξ::Float64
+end
+
+LogUtility() = LogUtility(1.0)
+
+(u::LogUtility)(c::Float64) =
+    ifelse(c > 1e-10, u.ξ*log(c), u.ξ*(log(1e-10) + 1e10*(c - 1e-10)))
+derivative(u::LogUtility, c::Float64) =
+    ifelse(c > 1e-10, u.ξ / c, u.ξ*1e10)
+
+"""
+Type used to evaluate CRRA utility. CRRA utility takes the form
+
+u(c) = ξ c^(1 - γ) / (1 - γ)
+
+Additionally, this code assumes that if c < 1e-10 then
+
+u(c) = ξ (1e-10^(1 - γ) / (1 - γ) + 1e-10^(-γ) * (c - 1e-10))
+"""
+struct CRRAUtility <: AbstractUtility
+    γ::Float64
+    ξ::Float64
+
+    function CRRAUtility(γ, ξ=1.0)
+        if abs(γ - 1.0) < 1e-8
+            error("Your value for γ is very close to 1... Consider using LogUtility")
+        end
+
+        return new(γ, ξ)
+    end
+end
+
+(u::CRRAUtility)(c::Float64) =
+    ifelse(c > 1e-10,
+           u.ξ * (c^(1.0 - u.γ) - 1.0) / (1.0 - u.γ),
+           u.ξ * ((1e-10^(1.0 - u.γ) - 1.0) / (1.0 - u.γ) + 1e-10^(-u.γ)*(c - 1e-10)))
+derivative(u::CRRAUtility, c::Float64) =
+    ifelse(c > 1e-10, u.ξ * c^(-u.γ), u.ξ*1e-10^(-u.γ))
+
+
+# Labor Utility
+
+"""
+Type used to evaluate constant Frisch elasticity (CFE) utility. CFE
+utility takes the form
+
+v(l) = ξ l^(1 + 1/ϕ) / (1 + 1/ϕ)
+
+Additionally, this code assumes that if l < 1e-10 then
+
+v(l) = ξ (1e-10^(1 + 1/ϕ) / (1 + 1/ϕ) - 1e-10^(1/ϕ) * (1e-10 - l))
+"""
+struct CFEUtility <: AbstractUtility
+    ϕ::Float64
+    ξ::Float64
+
+    function CFEUtility(ϕ, ξ=1.0)
+        if abs(ϕ - 1.0) < 1e-8
+            error("Your value for ϕ is very close to 1... Consider using LogUtility")
+        end
+
+        return new(ϕ, ξ)
+    end
+end
+
+(u::CFEUtility)(l::Float64) =
+    ifelse(l > 1e-10,
+           -u.ξ * l^(1.0 + 1.0/u.ϕ)/(1.0 + 1.0/u.ϕ),
+           -u.ξ * (1e-10^(1.0 + 1.0/u.ϕ)/(1.0 + 1.0/u.ϕ) + 1e-10^(1.0/u.ϕ) * (l - 1e-10)))
+derivative(u::CFEUtility, l::Float64) =
+    ifelse(l > 1e-10, -u.ξ * l^(1.0/u.ϕ), -u.ξ * 1e-10^(1.0/u.ϕ))
+
+
+"""
+Type used to evaluate elliptical utility function. Elliptical utility takes form
+
+v(l) = b (1 - l^μ)^(1 / μ)
+"""
+struct EllipticalUtility <: AbstractUtility
+    b::Float64
+    μ::Float64
+end
+
+# These defaults are pulled straight from Evans Phillips 2017
+EllipticalUtility(;b=0.5223, μ=2.2926) = EllipticalUtility(b, μ)
+
+(u::EllipticalUtility)(l::Float64) =
+    u.b * (1.0 - l^u.μ)^(1.0 / u.μ)
+derivative(u::EllipticalUtility, l::Float64) =
+    -u.b * (1.0 - l^u.μ)^(1.0/u.μ - 1.0) * l^(u.μ - 1.0)
+

test/runtests.jl

@@ -19,6 +19,7 @@ tests = [
         "markov_approx",
         "matrix_eqn",
         "mc_tools",
+        "modeltool",
         "quad",
         "quadsum",
         "random_mc",

test/test_ddp.jl

@@ -36,6 +36,8 @@ Tests for markov/ddp.jl
     Q_sa[3, :] = Q[2, 1, :]
     ddp0_sa = DiscreteDP(R_sa, Q_sa, beta, s_indices, a_indices)
 
+    @test issparse(ddp0_sa.Q)
+
     # List of ddp formulations
     ddp0_collection = (ddp0, ddp0_sa)
 

test/test_modeltool.jl

@@ -0,0 +1,65 @@
+@testset "Testing modeltools.jl" begin
+
+    @testset "Log Utility" begin
+    # Test log utility
+    u = LogUtility(2.0)
+    @test isapprox(u(1.0), 0.0)
+    @test isapprox(2.0*log(2.5), u(2.5))
+    @test isapprox(u.(ones(5)), zeros(5))
+    @test u(2.0) > u(1.0)  # Ensure it is increasing
+    @test isapprox(derivative(u, 1.0), 2.0)
+    @test isapprox(derivative(u, 2.0), 1.0)
+
+    # Test "extrapolation evaluations"
+    @test isapprox(u(0.5e-10), u(1e-10) + derivative(u, 1e-10)*(0.5e-10 - 1e-10))
+
+    @test isapprox(LogUtility().ξ, 1.0)  # Test default constructor
+
+    end
+
+    @testset "CRRA Utility" begin
+    # Test CRRA utility
+    u = CRRAUtility(2.0)
+    @test isapprox(u(1.0), 0.0)
+    @test isapprox((2.5^(-1.0) - 1.0) / (-1.0), u(2.5))
+    @test isapprox(u.(ones(5)), zeros(5))
+    @test u(5.0) > u(3.0)  # Ensure it is increasing
+    @test isapprox(derivative(u, 1.0), 1.0)
+    @test isapprox(derivative(u, 2.0), 0.25)
+
+    # Test "extrapolation evaluations"
+    @test isapprox(u(0.5e-10), u(1e-10) + derivative(u, 1e-10)*(0.5e-10 - 1e-10))
+
+    @test_throws ErrorException CRRAUtility(1.0)  # Test error throwing at γ=1.0
+    end
+
+    @testset "CFE Utility" begin
+    # Test CFE Utility
+    v = CFEUtility(2.0)
+    @test isapprox(v(1.0), -2.0/3.0)
+    @test isapprox(v(0.5), -v.ξ * 0.5^(1.0 + 1.0/v.ϕ) / (1.0 + 1.0/v.ϕ))
+    @test v(0.5) > v(0.85)  # Ensure it is decreasing
+    @test isapprox(derivative(v, 0.25), -0.5)
+    @test isapprox(derivative(v, 1.0), -1.0)
+
+    # Test "extrapolation evaluations"
+    @test isapprox(v(0.5e-10), v(1e-10) + derivative(v, 1e-10)*(0.5e-10 - 1e-10))
+
+    @test_throws ErrorException CRRAUtility(1.0)  # Test error throwing at ϕ=1.0
+    end
+
+    @testset "Elliptical Utility" begin
+    # Test Elliptical Utility
+    u = EllipticalUtility(1.0, 2.0)
+    @test isapprox(u(sqrt(0.75)), 0.5)
+    @test isapprox(u(sqrt(0.51)), 0.7)
+    @test u(0.5) > u(0.95)  # Ensure it is decreasing
+    @test isapprox(derivative(u, 0.0), 0.0)
+    @test isapprox(derivative(u, sqrt(0.5)), -1.0)
+
+    # Test default values
+    @test isapprox(EllipticalUtility().b, 0.5223)
+    @test isapprox(EllipticalUtility().μ, 2.2926)
+    end
+
+end