matlab - Find maximum of a utility function U(h1, h2)

Question

I have a matlab function that returns the utility value U(h1, h2) that depends on the number of work hours for two persons h1 and h2. Each person can work at maximum 3500 hours.

How do I write this code in the best way? For only one person I only loop the function from input value 0 to 3500 into a matrix and find the max value in the matrix. But now I have two values that can take the value 0 to 3500 and want to find the max output from the function based on these values.

I would be very grateful for any help!

edit: Sorry for late return! I haven't had time to look in to this until now. Anyway, here is the code:

function Test

global LocalTaxRate
global StateTax1Rate
global StateTax2Rate

LocalTaxRate = 0.3212; %Average 2017
StateTax1Rate = 0.25;
StateTax2Rate = 0.05;

h = [0 0]; % start value
lb = [0 0]; % lower bound of h
ub = [3500 3500]; % upper bound of h
myFun = @(h) -CalcUFamily(h(1), h(2)); % function to minimize with one input
Uoptimal = fmincon(myFun, [1000 1000], [], [], [], [], lb, ub)

end

function UFamily = CalcUFamily(hh,hw) %h = male, w = wife

(bunch of code until this step, which is the relevant one)

 YearlyWorkIncomeHusband = WageHH * hh;
    C = YearlyWorkIncomeHusband + MonthlyTaxableTransfers * 12 - CalcTaxSwe(YearlyWorkIncomeHusband + MonthlyTaxableTransfers * 12, YearlyWorkIncomeHusband, Age) + MonthlyNonTaxableTransfers * 12; %Netincome husband
    YearlyWorkIncomeWife = WageHW * hw;
    C = C + YearlyWorkIncomeWife + MonthlyTaxableTransfers * 12 - CalcTaxSwe(YearlyWorkIncomeWife + MonthlyTaxableTransfers * 12, YearlyWorkIncomeWife, Age) + MonthlyNonTaxableTransfers * 12; %Netincome husband + wife

    UFamily = alpha1 * log10(C/100000) + alpha11 * (log10(C/100000)^2) ...
        + alpha2 * (log10(T/1000-hh/1000)) + alpha22 * (log10(T/1000-hh/1000))^2 + alpha12 * log10(C/100000) * (log10(T/1000-hh/1000)) * 2 ... %husband
        + alpha3 * (log10(T/1000-hw/1000)) + alpha33 * (log10(T/1000-hw/1000))^2 + alpha13 * log10(C/100000) * (log10(T/1000-hw/1000)) * 2 ... %wife
        + alpha23 * (log10(T/1000-hh/1000)) * (log10(T/1000-hw/1000)) ... %common leisure parameter
        - alpha4 * Psa - bfc * Dw; %Dummies
end

Answer 1

You can use two approaches two solve your problem:

1. Generate a 2D mesh and find the global maximum

   [h1, h2] = meshgrid(0:100:3500);
   Uoptimal = max(max(U(h1, h2))
   

2. Use a multivariate minimisation solver, such as fmincon

   h = [1000 1000]; % start value
   lb = [0 0]; % lower bound of h
   ub = [3500 3500]; % upper bound of h
   Uoptimal = fmincon(@(h1, h2) -U(h1, h2), h, [], [], [], [], lb, ub);
   

   Note that I minimise -U, which is the same as maximising U.

Approach 1 will give you the global maxima, but is rather slow for obtaining an accurate result. Approach 2 is very accurate, but may get stuck at a local maxima, depending on the start value and the convexity of your utility function U.