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

**URL:** <https://discourse.julialang.org/t/solving-a-neoclassical-growth-model-via-pathsolver-but-solution-is-wrong/104979>\
**Category:** Optimization (Mathematical)\
**Created:** [October 15, 2023, 3:14pm UTC](https://discourse.julialang.org/t/solving-a-neoclassical-growth-model-via-pathsolver-but-solution-is-wrong/104979 "2023-10-15T15:14:40Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jovansam](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jovansam/32/15023_2.png) [@jovansam](https://discourse.julialang.org/u/jovansam)\
**Post date:** [October 15, 2023, 3:14pm UTC](https://discourse.julialang.org/t/solving-a-neoclassical-growth-model-via-pathsolver-but-solution-is-wrong/104979/1 "2023-10-15T15:14:40Z")

</div>

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)](https://sites.nd.edu/esims/files/2023/05/neoclassical_growth_notes_sp17.pdf)

```julia

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′]

```

---

<div class="post-metadata">

**Author:** ![jovansam](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jovansam/32/15023_2.png) [@jovansam](https://discourse.julialang.org/u/jovansam)\
**Post date:** [October 15, 2023, 3:44pm UTC](https://discourse.julialang.org/t/solving-a-neoclassical-growth-model-via-pathsolver-but-solution-is-wrong/104979/2 "2023-10-15T15:44:18Z")

</div>

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

```julia

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′]

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 15, 2023, 9:51pm UTC](https://discourse.julialang.org/t/solving-a-neoclassical-growth-model-via-pathsolver-but-solution-is-wrong/104979/3 "2023-10-15T21:51:23Z")

</div>

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