Suppose you have an Aiyagari model with endogenous labor supply. Under some restrictions on preferences, one might be able to solve for optimal labor conditional on (a’,a,z) in closed form. Then one could easily write the return function with (a’,a,z) as the only input arguments (besides scalar parameters) and obtain n using the closed form solution.
However, for many popular utility functions, the closed form solution is not available and one has to solve numerically the first order condition. It turns out that you cannot call fzero or fsolve inside the return function (due to limitations in function calls by arrayfun). However, you (or your AI) can easily write a rootfinding loop such as bisection and viola’!
function [F, n] = static_problem(ap, a, z, ...
w, r, sigma, phi, nu, nmax)
% Resources excluding labor income
x = (1+r)*a - ap;
wage = w*z;
% Ensure positive consumption
cmin = 1e-10;
% Minimum feasible labor
nlo = max(0.0, (cmin - x)/wage);
nhi = nmax;
% Infeasible even at maximum labor
if nlo >= nmax
n = nmax;
F = -Inf;
return
end
% FOC at lower bound
c = x + wage*nlo;
glo = wage*c^(-sigma) - phi*nlo^nu;
% FOC at upper bound
c = x + wage*nhi;
ghi = wage*c^(-sigma) - phi*nhi^nu;
% Lower corner
if glo <= 0
n = nlo;
% Upper corner
elseif ghi >= 0
n = nhi;
else
% Interior solution: bisection
for iter = 1:40
nmid = 0.5*(nlo+nhi);
c = x + wage*nmid;
gmid = wage*c^(-sigma) ...
- phi*nmid^nu;
if gmid > 0
nlo = nmid;
else
nhi = nmid;
end
end
n = 0.5*(nlo+nhi);
end
% Consumption at optimal labor
c = x + wage*n;
% Current utility
if sigma == 1
uc = log(c);
else
uc = (c^(1-sigma)-1)/(1-sigma);
end
un = phi*n^(1+nu)/(1+nu);
F = uc - un;
end
This is faster and more precise that computing optimal labor on a discrete grid. The user can just define a model like this as a model with exogenous labor (i.e. without d variable)
