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

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

``````

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, 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

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