Any optimization libraries/guides?

Good evening,

First of all, I know this problem is unrelated to the lectures but this is the first forum that comes to mind when thinking about Python and economics.

I’m currently trying to solve an optimization problem using Pyomo. Specifically, I’m trying to solve a maximization problem subject to constraints but I’m having trouble understanding the structure of Pyomo’s code. Also, I cannot seem to find any relevant guides or examples applied to economics.

Are there any other optimization libraries for Python with such guides?

I’m not familiar with Julia but I’ve seen JuMP allows mathematical optimization. Any resources on where to start?

Thank you for your time.

Kind regards,
Tomás.

1 Like

Hi @ta.rodriguez,

Good question. Can you tell us a bit more about your problem. Is it convex, for example?

I’m a bit out of touch because I typically use dynamic programming with a low dimensional action space and rely on scipy solvers.

@cc7768 , do you have any thoughts on this?

1 Like

Hi @john.stachurski ,

You’re always helpful.

As I’ve said in a previous post. I’m currently working on my master’s thesis which is an analysis of the rural credit markets in developing economies.

I have a single agent problem where a representative entrepreneur seeks to maximize her life-time utility under the following CRRA function:

\max \sum_{t=0}^{\infty}{\beta^{t}(\frac{C_{t}^{1-\theta}-1}{{1-\theta}})}

She has access to a decreasing returns technology that only uses capital K to produce a consumption good Y:

g(K_{t})=zK_{t}^{\alpha}

Where z is a permanent ability component and \alpha<1 is the degree of decreasing returns to capital.

At time t, the entrepreneur must decide how much capital she wants to have in the next period K_{t+1}.

She must also decide how much credit she wants to take, which can be a combination of formal credit F_{t+1}, and/or informal credit I_{t+1}.

If she takes no credit, then she will self-finance her production.

The formal credit market imposes a credit constraint on the entrepreneur:

qK_{t}\geq F_{t+1}

This constraint impedes the formal lender to grant more than a proportion q\in [0,1] of the capital owned by the entrepreneur at time t.

The informal credit market does not impose any constraints on the entrepreneur.

Thus, the maximization problem is the following:

\max_{[C_{t}, K_{t+1}, I_{t+1}, F_{t+1}]} \sum_{t=0}^{\infty}{\beta^{t}(u(C_{t}))}
s.t.
g(K_{t})+(1-\delta)K_{t}+I_{t+1}+F_{t+1} \geq C_{t}+K_{t+1}+(1+i_{t})I_{t}+(1+f_{t})F_{t}
K_{t}\geq 0
I_{t+1}\geq 0
F_{t+1}\geq 0
qK_{t}\geq F_{t+1}

And the Lagrangian function:

L=\sum_{t=0}^{\infty}{\beta^{t}}\{u(C_{t})+\lambda_{t}[g(K_{t})+(1-\delta)K_{t}
+I_{t+1}+F_{t+1}-C_{t}-K_{t+1}-(1+i_{t})I_{t}+(1+f_{t})F_{t}]
+\mu_{t}(K_{t})+\gamma_{t}(I_{t+1})+\sigma_{t}(F_{t+1})+\psi_{t}(qK_{t}-F_{t+1})\}

This is a concave problem since u(C_{t}) and g(K_{t}) are both concave and the rest of the arguments are linear. On a previous post, I was trying to find the steady state of this problem. My advisor suggested that I should focus on the solution of the problem instead.

Again, thank you for your help and time.

1 Like

I’m in the same boat as @john.stachurski, my optimization problems are typically relatively low dimensional. I might consider trying to discretize the choice variables first if you only have 3-4 of them – This should give you a relatively accurate solution (and has the benefit of being “bulletproof” in spite of being a bit slow).

I’ve also been meaning to experiment more with Bayesian optimization (see, this package) for some of optimization problems I’ve worked on, but I haven’t gotten around to it yet.

1 Like
On a previous post, I was trying to find the steady state of this problem. My advisor suggested that I should focus on the solution of the problem instead.

Your advisor is correct. The SS is trivial as the constraints won’t be binding (well, qK\geqslant F might be binding). The whole dynamics is a bit more interesting than that because it’s an “occasional binding constraint”.
To solve this you can follow a similar approach as in Aiyagari model.

1 Like

@ta.rodriguez I was looking for similar info.
I found this lecture by Hua Zhou very helpful.
First the flowchart:


Next a table w/ some solvers in R/Matlab/Julia/Python:

One package I find useful that is not on this list in Optim.jl.

I think it would be great for QuantEcon lectures to include at least one lecture on the optimization ecosystem…

Thanks @azev77, that’s super useful.

I think it would be great for QuantEcon lectures to include at least one lecture on the optimization ecosystem…

I fully agree but I’m tied up with other priorities right now :grimacing: (R&Rs plus writing a hardcopy textbook on macroeconomic modeling).

1 Like

To be clear, there already is a lecture on optimization solver packages, my suggestion is to incorporate the structure in the flowchart above into the lecture: (is it convex? then use Convex.jl etc …)

@ta.rodriguez
You have a neoclassical growth model (NCGM) w/ two types of debt:

  1. Collateralized debt: F_{t} at interest rate f_{t} s.t. constraint F_{t} \leq q K_{t}
  2. Unsecured debt: I_{t} at interest rate i_{t}
    We need i_{t} > f_{t} for this problem to be interesting.
    It makes sense: unsecured credit card rates (20%) are higher than mortgage rates (3%)

Without debt you have a routine NCGM w/ a known closed form solution (Brock-Mirman) w/ u(c)=log(c), \delta=1, g(k)=z k^\alpha.
This should be your starting point.

The easiest way I know how to solve this problem is to cast it in continuous time & use an Optimal Control package. I use NLOptControl.jl which can handle inequality constraints, but currently only allows finite horizon.

J = \max \int_{t=0}^{t=T} e^{-\rho t} u(c_{t}) dt
s.t.
\dot{k}_{t} = k_{t}^{\alpha} -\delta k_{t} - c_{t} - pmt_{1t} - pmt_{2t}
\dot{b}_{1t} = r_{1} \times b_{1t} - pmt_{1t}
\dot{b}_{2t} = r_{2} \times b_{2t} - pmt_{2t}
b_{1t} \leq q \times k_{t}
c_{t}, k_{t}, b_{1t}, b_{2t} \geq 0
k_{0} given, k_{T}=0
b_{i0}=0, b_{iT}=0 for i\in \{1,2\} (born w/o debt, die w/o debt)

State variables: k_t, b_{1t}, b_{2t}
Choice variables: c_t, pmt_{1t}, pmt_{2t}
Parameters: T, \rho, u(c), \alpha, \delta, r_{1}, r_{2}, q

Case 1: no borrowing

using NLOptControl

# NCGM No Borrowing
n=define(
    numStates=1, numControls=1,
    X0=[5.0], XF=[0.01],
    XL=[0.0], XU=[NaN],   # k(t) >= 0.0
    CL=[0.01], CU=[NaN]   # c(t) >= 0.0
    )
states!(n,[:k];descriptions=["k(t)"])
controls!(n,[:c];descriptions=["c(t)"])
dx=[
    :( (k[j]^(.67)) - 0.3 * k[j] - c[j] )
    ]
dynamics!(n,dx)
configure!(n;
    (:integrationScheme=>:trapezoidal),
    (:finalTimeDV=>false),(:tf=>100.0))
# @NLconstraint
obj=integrate!(n, :( exp(-.01 * j) * log(c[j]) ) );
@NLobjective(n.ocp.mdl, Max, obj)
optimize!(n)
allPlots(n)

Case 2: Collateralized Borrowing

# NCGM + Collateralized Borrowing
n=define(
    numStates=2, numControls=2,
    X0=[5.0,0.0], XF=[0.01,0.0],
    XL=[0.0,NaN], XU=[NaN,NaN],   # k(t) >= 0.0
    CL=[0.01,-20.0], CU=[NaN,NaN]   # c(t) >= 0.0
    )
states!(n,[:k,:b1];descriptions=["k(t)","b1(t)"])
controls!(n,[:c,:pmt1];descriptions=["c(t)","pmt1(t)"])
dx=[
    :( (k[j]^(.67)) - 0.3 * k[j] - c[j] - pmt1[j]),
    :( 0.01 * b1[j] - pmt1[j] )
    ]
dynamics!(n,dx)
configure!(n;
    (:integrationScheme=>:trapezoidal),
    (:finalTimeDV=>false),(:tf=>100.0)
    )
# @NLconstraint
obs_con = @NLconstraint(n.ocp.mdl,
    [i=1:n.ocp.state.pts-1],
    n.r.ocp.x[i,2] <= 0.25 * n.r.ocp.x[i,1]  #b[t] <= 0.5*k[t]
    )
newConstraint!(n,obs_con,:obs_con)
#
obj=integrate!(n, :( exp(-.01 * j) * log(c[j]) ) );
@NLobjective(n.ocp.mdl, Max, obj)
optimize!(n)
allPlots(n)

Case 3: Collateralized Borrowing + Unsecured debt

# NCGM + Collateralized Borrowing + Unsecured debt
n=define(
    numStates=3, numControls=3,
    X0=[5.0,0.0,0.0], XF=[0.01,0.0,0.0],
    XL=[0.0,NaN,0.0], XU=[NaN,NaN,NaN],   # k(t) >= 0.0, no saving
    CL=[0.01,-20.0,-20.0], CU=[NaN,NaN,NaN]   # c(t) >= 0.0
    )
states!(n,[:k,:b1,:b2];descriptions=["k(t)","b1(t)","b2(t)"])
controls!(n,[:c,:pmt1,:pmt2];descriptions=["c(t)","pmt1(t)","pmt2(t)"])
dx=[
    :( (k[j]^(.67)) - 0.3 * k[j] - c[j] - pmt1[j] - pmt2[j]),
    :( 0.01 * b1[j] - pmt1[j] ),
    :( 0.03 * b2[j] - pmt2[j] )
    ]
dynamics!(n,dx)
configure!(n;
    (:integrationScheme=>:trapezoidal),
    (:finalTimeDV=>false),(:tf=>100.0)
    )
# @NLconstraint
obs_con = @NLconstraint(n.ocp.mdl,
    [i=1:n.ocp.state.pts-1],
    n.r.ocp.x[i,2] <= 0.25 * n.r.ocp.x[i,1]  #b[t] <= 0.5*k[t]
    )
newConstraint!(n,obs_con,:obs_con)
#
obj=integrate!(n, :( exp(-.01 * j) * log(c[j]) ) );
@NLobjective(n.ocp.mdl, Max, obj)
optimize!(n)
allPlots(n)

You’re gonna have to better tune the solver parameters for more precise solutions…
Good luck

2 Likes

@azev77,

Thank you for all your support!

My response is a bit late because I had to deliver my final document to my faculty. I did not use your code and ended up finding policy functions in MatLab. I can share the document and/or code if anyone is interested.

Before I did this I was able to develop a code to find optimal values for this problem when you give the code some information, which I want to share if somebody ever wants to solve a similar problem.

#!conda install -c conda-forge pyomo -y
#!conda install -c conda-forge ipopt -y

import pyomo
import pyomo.opt
import pyomo.environ as pyo
from pyomo.opt import SolverFactory
from pyomo.environ import *

import pandas

from numba import njit, jitclass, float64
import numpy as np
import matplotlib.pyplot as plt

data = {
    ('1111'): {'C': 50, 'K': 15, 'F': 0, 'f': 0, 'F_next': 0, 'I': 0, 'I_next': 0, 'K_next': 25, 'i': 0},
    ('1112'): {'C': 100, 'K': 60, 'F': 20, 'f': 0.05, 'F_next': 22, 'I': 0, 'I_next': 0, 'K_next': 75, 'i': 0},
    ('1112'): {'C': 100, 'K': 60, 'F': 18, 'f': 0.07, 'F_next': 22, 'I': 0, 'I_next': 0, 'K_next': 75, 'i': 0},
    ('1113'): {'C': 110, 'K': 65, 'F': 22, 'f': 0.06, 'F_next': 26, 'I': 0, 'I_next': 0, 'K_next': 75, 'i': 0},
    ('1114'): {'C': 60, 'K': 20, 'F': 0, 'f': 0, 'F_next': 0, 'I': 10, 'I_next': 20, 'K_next': 30, 'i': 0.25},
    ('1115'): {'C': 90, 'K': 55, 'F': 12, 'f': 0.06, 'F_next': 15, 'I': 8, 'I_next': 150, 'K_next': 65, 'i': 0.20},
}

θ = 2
β = 0.958
δ = 0.075
α = 0.85
z = 2
q = 0.5

Obs = data.keys()

for s in Obs:

    model = pyo.ConcreteModel()
    model.K_next = pyo.Var(Obs, within=pyo.NonNegativeReals)
    model.I_next = pyo.Var(Obs, within=pyo.NonNegativeReals)
    model.F_next = pyo.Var(Obs, within=pyo.NonNegativeReals)
    
    model.U = pyo.Objective(expr = log(z * data[s]['K'] ** α + (1 - δ) * data[s]['K'] + model.I_next[s] + model.F_next[s] - model.K_next[s] - (1 + data[s]['i']) * data[s]['I'] - (1 + data[s]['f']) * data[s]['F']), sense=pyo.maximize)

    def Capital_rule(model):
        return (0, data[s]['K'] >= 0, None)
    model.Capital = Constraint(rule = Capital_rule)

    def Informal_rule(model):
        return (0, data[s]['I_next'] >= 0, None)
    model.Informal = Constraint(rule = Informal_rule)

    def Formal_rule(model):
        return (0, data[s]['F_next'] >= 0, None)
    model.Formal = Constraint(rule = Formal_rule)

    def Collateral_rule(model):
        return (0, q * data[s]['K'] >= data[s]['F_next'], None)
    model.Collateral = Constraint(rule = Collateral_rule)

    solver = pyo.SolverFactory('ipopt')
    solver.options['max_iter'] = 1000
    solver.solve(model)
    print()
    print('Observation = ', s)
    print()
    print('Utility = ', pyo.value(model.U))
    print()
    print('Capital = ', pyo.value(model.Capital))
    print()
    print('Informal = ', pyo.value(model.Informal))
    print()
    print('Formal = ', pyo.value(model.Formal))
    print()
    print('Collateral = ', pyo.value(model.Collateral))
    print()
    print('K_next = ', pyo.value(model.K_next[s]))
    print()
    print('F_next = ', pyo.value(model.F_next[s]))
    print()
    print('I_next = ', pyo.value(model.I_next[s]))
    print()
    
    data_estrella[s] = {
    (s): {'Utility': pyo.value(model.U), 'Capital': pyo.value(model.Capital), 'Informal': pyo.value(model.Informal), 'Formal': pyo.value(model.Formal), 'Collateral': pyo.value(model.Collateral), 'K_next': pyo.value(model.K_next[s]), 'I_next': pyo.value(model.I_next[s]), 'F_next': pyo.value(model.F_next[s])}
    }

print(data_estrella)