Can't seem to get NLsolve to converge

New to Julia
1

Please note: I’ve cross posted this question on StackOverflow.
here

I’m trying to solve a life cycle problem in economics using Julia but I’m having trouble with NLsolve. The model boils down to trying to solve two a two equation system to find optimal leisure hours and capital stock each working period. The economic agent after retirement sets leisure = 1 and I only need to solve a single non linear equation for capital. This part works fine. It’s solving the two equation system that seems to break down.

As I’m fairly new to Julia / programming in general so any advice would be very helpful. Also advice / points / recommendations on all aspects of the code will be greatly appreciated. The model is solved backwards from the final time period.

My attempt

using Parameters
using Roots
using Plots
using NLsolve
using ForwardDiff

Model = @with_kw (α = 0.66,
                  δ = 0.02,
                  τ = 0.015,
                  β = 1/1.01,
                  T = 70,
                  Ret = 40,
                  );

function du_c(c, l, η=2, γ=2)
    if c>0 && l>0
        return (c+1e-6)^(-η) * l^((1-η)*γ)
    else
        return Inf
    end
end

function du_l(c, l, η=2, γ=2)
    if l>0 && c>0
        return γ * (c+1e-6)^(1-η) * l^(γ*(1-η)-1)
    else
        return Inf
    end
end

function create_euler_work(x, y, m, k, l, r, w, t)
    # x = todays capital, y = leisure
    @unpack α, β, τ, δ, T, Ret = m
    c_1 = x*(1+r) + (1-τ)*w*(1-y) - k[t+1]
    c_2 = k[t+1]*(1+r) + (1-τ)*w*(1-l[t+1]) - k[t+2]
    return du_c(c_1,y) - β*(1+r)*du_c(c_2,l[t+1])
end

function create_euler_retire(x, m, k, r, b, t)
    # Holds at time periods Ret onwards 
    @unpack α, β, τ, δ, T, Ret = m
    c_1 = x*(1+r) + b - k[t+1]
    c_2 = k[t+1]*(1+r) + b - k[t+2]
    return du_c(c_1,1) - β*(1+r)*du_c(c_2,1)
end

function create_euler_lyw(x, y, m, k, r, w, b, t)
    # x = todays capital, y = leisure
    @unpack α, β, τ, δ, T, Ret = m
    c_1 = x*(1+r) + (1-τ)*w*(1-y) - k[t+1]
    c_2 = k[t+1]*(1+r) + b - k[t+2]
    return du_c(c_1,y) - β*(1+r)*du_c(c_2,1)
end


function create_foc(x, y, m, k, r, w, t)
    # x = todays capital, l= leisure
    @unpack α, β, τ, δ, T = m
    c = x*(1+r) + (1-τ)*w*(1-y) - k[t+1]
    return du_l(c,y) - (1-τ)*w*du_c(c,y)
end


function life_cycle(m, guess, r, w, b, initial)
    @unpack α, β, τ, δ, T, Ret = m
    k = zeros(T+1);
    l = zeros(T);
    k[T] = guess

    println("Period t = $(T+1) Retirment, k = $(k[T+1]), l.0 = NA")
    println("Period t = $T Retirment, k = $(k[T]), l = 1.0")
    ########################## Retirment ################################

    for t in T-1:-1:Ret+1 
        euler(x) = create_euler_retire(x, m, k, r, b, t)
        k[t] = find_zero(euler, (0,100))
        l[t] = 1
        println("Period t = $t Retirment, k = $(k[t]), l = $(l[t])")
    end 
    ###################### Retirement Year  #############################
    for t in Ret:Ret
        euler(x,y) = create_euler_lyw(x, y, m, k, r, w, b, t)
        foc(x,y) = create_foc(x, y, m, k, r, w, t)

        function f!(F, x)
            F[1] = euler(x[1], x[2])
            F[2] = foc(x[1], x[2])
        end

        res = nlsolve(f!, [5; 0.7], autodiff = :forward)
        k[t] = res.zero[1]
        l[t] = res.zero[2]
        println("Period t = $t Working, k = $(k[t]), l = $(l[t])")
    end 
    ############################ Working ###############################
    for t in Ret-1:-1:1
        euler(x,y) = create_euler_work(x, y, m, k, l, r, w, t)
        foc(x,y) = create_foc(x, y, m, k, r, w, t)

        function f!(F, x)
            F[1] = euler(x[1], x[2])
            F[2] = foc(x[1], x[2])
        end

        res = nlsolve(f!, [5; 0.7], autodiff = :forward)
        k[t] = res.zero[1]
        l[t] = res.zero[2]
        println("Period t = $t Working, k = $(k[t]), l = $(l[t])")
    end
    #####################################################################
    return k[1] - initial, k, l
end


m = Model();

residual, k, l = life_cycle(m, 0.3, 0.03, 1.0, 0.0, 0.0)

The code seems to break on period 35 with the error “During the resolution of the nonlinear system, the evaluation of following equations resulted in a non-finite number: [1,2]” However the solutions seem to go weird at period 37.

2 Likes
2

I think you are running into non-positive consumption or labor choices, get infinities, which then combine into a NaN.

I would recommend reformulating the problem in a way that either avoids or incorporates this, eg using a reparametrization \tilde{l} = \log(l) and \tilde{c} = \log(c) and solving for all periods simultaneously (as a nonlinear system). The infinite-lifetime (dynasty) model or a one-shot version should provide a reasonable initial guess for the solver.

2 Likes
3

Thanks for the reply. Your exactly right! Turns out that instead of declaring c,l<0 = Inf I can just set it to some arbitrary large number 1e6 or so. However it only converges for a very narrow interval of guesses. Not sure why that is. Cheers.

4

I didn’t want to formulate the problem as a large non linear system. With the way it’s written it exploits the time structure to make it more efficient, well hopefully.

5

I don’t think so; breaking it up to T intertemporal problems should be less efficient IMO because the solver will be constrained on the manifold defined by the Euler equations.

6

I’ll certainly try your suggestion and report the results. Will be an Interesting comparison. I’m relatively new to programming in general so all advice is appreciated. Thanks.