Solving a neoclassical growth model via PATHSolver (but solution is wrong)

Specific Domains Optimization (Mathematical)
1

Hi. I am solving the following neoclassical growth model via the PATHSolver (i.e. MCP formulation). I know the solution is wrong as consumption is unusually high and not inline with what I get from GAMS. I suspect this stems from how I am specifying the complementarity condition. I would appreciate any guidance to have the correct formulation (and hence an answer that makes sense.) This type of model can be found here neoclassical_growth_notes_sp17.pdf (nd.edu)


print("\033c")
print("#------------------------------------\n")
print("country-level SCC\n")
print("#------------------------------------\n")

using JuMP
using PATHSolver

function scc()
    η = 2
    ρ = 0.001
    β = 1/(1+ρ)
    T = 10 # number of periods
    t = 1:T
    α = 0.3
    δ_k = 0.1
    K₀ = 5.0

    L = ones(size(t))
    Aₖ = ones(size(t))

    u(c,L) = L*(c/L^(1-η)) / (1-η) - 1
    ∂u∂c(c,L) = (c/L)^(-η)

    scc = Model(PATHSolver.Optimizer)
    
    @variable(scc, c[t]>=0, start = 0.3)
    @variable(scc, K[t]>=0, start = K₀)
    @variable(scc, λ[t]>=0, start = 1.0)
    
    fix(K[1], K₀, force=true)
    
    #@constraint(scc, [i=t], (c[i]/L[i])^(-η) - λ[i] ⟂ c[i])
    @constraint(scc, [i=t], -1*(∂u∂c(c[i],L[i]) - λ[i]) ⟂ c[i])
    @constraint(scc, [i=1:(T)], -1*(Aₖ[i]*K[i]^α*L[i]^(1-α) - c[i] + (1-δ_k)*K[i] - (i<T ? K[i+1] : K[i])) ⟂ λ[i])
    @constraint(scc, [i=1:(T-1)], -1*(-λ[i] + β*λ[i+1]*(Aₖ[i+1]*α*K[i+1]^(α-1)*L[i+1]^(1-α) + (1-δ_k))) ⟂ K[i+1])
    
    optimize!(scc)
    print(solution_summary(scc))

    y = zeros(T)
    y = Aₖ.*value.(K).^α.*L.^(1-α)
    k′ = y - value.(c) + (1-δ_k)*value.(k)
    return value.(c), value.(K), y, value.(λ), k′
end

c, k, y, λ, k′ = scc()

[c k y λ k′]
2

Solved it myself. Here is the correct code (if it helps anyone else)


print("\033c")
print("#------------------------------------\n")
print("country-level SCC\n")
print("#------------------------------------\n")

using JuMP
using PATHSolver

function scc()
    η = 2
    ρ = 0.001
    β = 1/(1+ρ)
    T = 10 # number of periods
    t = 1:T
    α = 0.3
    δ_k = 0.1
    K₀ = 5.0

    L = ones(size(t))
    Aₖ = ones(size(t))

    u(c,L) = L*(c/L^(1-η)) / (1-η) - 1
    ∂u∂c(c,L) = (c/L)^(-η)

    scc = Model(PATHSolver.Optimizer)
    
    @variable(scc, c[t]>=0, start = 0.3)
    @variable(scc, K[t]>=0, start = K₀)
    
    # formulation with shadow price on capital as free variable
    #@variable(scc, λ[t], start = 1.0)
    #@constraint(scc, [i=1:(T)], -1*(Aₖ[i]*K[i]^α*L[i]^(1-α) + (1-δ_k)*K[i] - c[i] - (i<T ? K[i+1] : K[i])) ⟂ λ[i])
    
    #formulation with shadow price on capital as a signed (i.e. positive) price
    @variable(scc, λ[t]>=0, start = 1.0)
    @constraint(scc, [i=1:(T)], (Aₖ[i]*K[i]^α*L[i]^(1-α) + (1-δ_k)*K[i] - c[i] - (i<T ? K[i+1] : K[i])) ⟂ λ[i])
    
    fix(K[1], K₀, force=true)
    
    #@constraint(scc, [i=t], (c[i]/L[i])^(-η) - λ[i] ⟂ c[i])
    @constraint(scc, [i=t], -1*(∂u∂c(c[i],L[i]) - λ[i]) ⟂ c[i])
    @constraint(scc, [i=1:(T-1)], -1*(-λ[i] + β*λ[i+1]*(Aₖ[i+1]*α*K[i+1]^(α-1)*L[i+1]^(1-α) + (1-δ_k))) ⟂ K[i+1])
    
    optimize!(scc)
    print(solution_summary(scc))

    y = zeros(T)
    y = Aₖ.*value.(K).^α.*L.^(1-α)
    k′ = y - value.(c) + (1-δ_k)*value.(k)
    return value.(c), value.(K), y, value.(λ), k′
end

c, k, y, λ, k′ = scc()

[c k y λ k′]
2 Likes
3

Glad you figured it out. It can be very tricky to debug models like this if you forget a single minus sign :smile:

1 Like