1

我正在尝试使用 Julia 解决经济学中的生命周期问题，但我在使用 NLsolve 时遇到了麻烦。该模型归结为尝试解决两个二等式系统，以找到每个工作期间的最佳休闲时间和资本存量。退休后的经济代理人设置休闲= 1，我只需要求解一个资本非线性方程。这部分工作正常。它正在解决似乎崩溃的两个方程系统。

因为我对 Julia / 编程很陌生，所以任何建议都会非常有帮助。此外，我们将不胜感激有关代码各个方面的建议/要点/建议。该模型从最后一个时间段向后求解。

我的尝试

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)

代码似乎在第 35 期中断，并出现错误“在非线性系统的解析过程中，对以下方程的评估导致了一个非有限数：[1,2]”但是在第 37 期，解决方案似乎变得很奇怪。

4

0 回答 0