Max a function wrt two variables

I have been asking this at another forum, and got nice help - but I still can't get it quite to work.

I have a function for calculating the utility of for working a certain number of hours (yearly) for a husband and wife in a family. I want to maximize the utility function, CalcUFamily, with respect to the number of working hours (1 to 3500) for the wife and husband, separately and give me the optimal work hours for the husband and wife.

I tried the code below, but can't get it to work properly. Can you help me?

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 irrelevant code skipped)
else
  Psa = 0;
  Dw = 1;
  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
end

Answer 1

[Image Analyst]

The code crashes at this point:

function UFamily = CalcUFamily(hh,hw) %h = male, w = wife
(bunch of irrelevant code skipped)
else

You need to start the else block off with an "if" statement. And the "(bunch..." needs to be in a comment - a line starting with %. So, try something like this:

function UFamily = CalcUFamily(hh,hw) %h = male, w = wife
% (bunch of irrelevant code skipped)
if someCondition
  % whatever...
else

[KGB91]

I do, I just skipped it out because it wasn't relevant for the code here, I thought. I have a long code of calculating the alpha-values and an if before the else. Do you want me to post more of the code?

[Image Analyst]

Well, that pretty much stopped me in my tracks and I quit after that. I'll leave it up to you.

[KGB91]

Let's make it easy, this simplification should do it:

function UFamily = CalcUFamily(hh,hw) %h = male, w = wife
    Psa = 0;
    Dw = 1;
    C = 100 * hh + 100 * hw;
T = 4000
Psa  = 0;
bfc = 0;
alpha1 = 8.893000000000000;
alpha11 = 2.257000000000000;
alpha12 = -0.42;
alpha13 = -0.630000000000000;
alpha2 = 18.554452000000005;
alpha22 = -12.594000000000000;
alpha23 = 0.857000000000000;
alpha3 = 3.161244000000003;
alpha33 = -12.503000000000000;
alpha4 = 12.469726000000000;
    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

[Walter Roberson]

What difficulty are you observing?

When I run your code, I get an exit flag of 1 from the fmincon, which is the exit flag that indicates success.

[Walter Roberson]

Note that you coded 0 as your lower bounds, but your problem description gives 1 as the lower bounds. The difference is important in this case, as your function has a maxima at points arbitrarily close to 0 (and is undefined at (0,0))

[KGB91]

The trouble is that I cannot get proper values for the optimal working time, which should be around 1240 for the husband and 2950 for the wife in this case. The code just gives me

A =

   1.0e-09 *

    0.0000    0.1782

B =

 -265.3649

if I replace Uoptimal with [A,B] in order to get the optimal working hours. I made a mistake in the problem description, it should go from 0 to 3500. But the 0 part of the code is not finished, so let's skip the 0 for now - I will add it later when the rest of the code works.

Answer 2

[Walter Roberson]

The maximum is at x1 = 1299.789403251401341605335676596 and x2 = 2993.8956914836407431811398740633 giving a maximum of 11.969999104694875886276737089086, compared to only 11.96200980280465830244970897474 for [1240, 2950]

syms x1 x2
F = CalcUFamily(x1,x2);
dx1 = diff(F,x1);
dx2 = diff(F,x2);
sol = vpasolve(dx1, dx2);
best_x1 = sol.x1; best_x2 = sol.x2;

However, there is a region close to 0 where in theory the function goes infinity, in pretty much a 1/x relationship (that is, the value at 10^-(N) is about 1/10th of the value at 10^-(N+1) ). It does not grow especially quickly -- at 10^-100000000 the value has only reached 22570000455447658.275295509433155 . None the less, it goes to infinity.

[KGB91]

Thanks :) However I cant get the code to work, it sets hh to x1 and hw to x2?

I know the issue with the function at hh = hw = 0, I am working on it. The parameters aren't finished yet.

[KGB91]

Sorry for asking again, but could you tell me how you got the solution? :)

[Walter Roberson]

Do you have the symbolic toolbox? If you do then the code is what I posted.

[KGB91]

I don't since I have Matlab R2014B. Thanks, I'll upgrade now and post the result here. :)

[Walter Roberson]

You could have the Symbolic toolbox in MATLAB R2014b; the latest version is not necessary, but installing the (licensed) toolbox would be.

[KGB91]

Yeah, it took longer to update than expected so that was probably something I should have thought about...

[KGB91]

I tried to update it to Matlab 2017a. The problem now it that function Test thinks that hh = x1 and hw = x2, which it cannot test (the wage is based on the number of working hours). The function test is now the following code (copied from you):

function Test

global LocalTaxRate
global StateTax1Rate
global StateTax2Rate

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

syms x1 x2
F = CalcUFamily(x1,x2);
dx1 = diff(F,x1);
dx2 = diff(F,x2);
sol = vpasolve(dx1, dx2);
best_x1 = sol.x1; best_x2 = sol.x2;

end

[Walter Roberson]

The LocalTaxRate, StateTax1Rate, and StateTax2Rate that you assign there are not used in your CalcUFamily code. They are irrelevant. So you simply call

syms x1 x2
F = CalcUFamily(x1,x2);
dx1 = diff(F,x1);
dx2 = diff(F,x2);
sol = vpasolve(dx1, dx2);
best_x1 = sol.x1; best_x2 = sol.x2;

None of the code you show needs to test anything.

If it will make it easier for you to understand, then use

syms hh hw
F = CalcUFamily(hh,hw);
dhh = diff(F,hh);
dhw = diff(F,hw);
sol = vpasolve(dhh, dhw);
best_hh = sol.hh; best_hw = sol.hw;

[KGB91]

Thanks, now I see the problem! I tried to shorten and simplify the code here i order to save some space. I shouldn't have done that, that is why it isn't working. Here is the real code:

function uFamily = CalcUFamily(hh,hw) %h = male, w = wife
BigCity = 1; %Always 1, does not exist in the data...
T = 4000; %Total amount of time
MaxWorkTime = 3500;

%These will later be input from micro data
HighestEducationStatus = 1;
Citizenship = 1;
Age = 22;
WageHH = 100; %husband
WageHW = 100; %wife
MonthlyTaxableTransfers = 0;
MonthlyNonTaxableTransfers = 0;
NumberOfChildren = 0;
MonthlyHousingCostTotal = 1000;

if HighestEducationStatus <= 2
    PrimarySchool = 1;
    SecondarySchool = 0;
else
    PrimarySchool = 0;
    SecondarySchool = 1;
end

if Citizenship == 1
    BornInTheCountry = 1; %Citizenship for Sweden, always 1 in the UK
else
    BornInTheCountry = 0;
end

if (Age > 17) && (Age < 30)
    Age18to29 = 1;
    Age30to34 = 0;
elseif (Age > 29) && (Age < 35)
  Age30to34 = 1;
  Age18to29 = 0;
else
  Age30to34 = 0;
  Age18to29 = 0;
end

alpha1 = 8.893;
alpha11 = 2.257;
alpha22 = -12.594;
alpha33 = -12.503;
alpha12 = -0.420;
alpha13 = -0.630;
alpha23 = 0.857;

%-----------COMMON---------------

%Welfare
SA1 = 21.970;
SA2 = -9.399;
SA3 = -12.571;

%---------Husband, alfa2-----------

%LEISURE
InterceptLeisure = 34.277;
BigCityLeisure = 0.369;
PrimarySchoolLeisure = -0.794;
SecondarySchoolLeisure = -1.235;
AgeEffectLeisure = -0.296;
AgeEffect2Leisure = 0.332;    

%FIXED COST
bfch = 10.542;

%WELFARE - alpha4
InterceptWelfare = 13.826;
BigCityWelfare = -0.094;
SwedishBornWelfare = 2.890;
PrimarySchoolWelfare = -1.134;
AgeEffectWelfare18to29 = 6.005;
AgeEffectWelfare30to34 = 2.536;

%UNOBSERVED HETEROGENEITY
L1 = 7.198; %Leisure
L2 = -10.940;
L3 = 3.742;

L = L1 * 0.014 + L2 * 0.966 + L3 * 0.02; %Parameter * probability
SA = SA1 * 0.014 + SA2 * 0.966 + SA3 * 0.02; %Parameter * probability
alpha2 = InterceptLeisure + BigCityLeisure * BigCity + PrimarySchoolLeisure * PrimarySchool + SecondarySchoolLeisure * SecondarySchool + AgeEffectLeisure * Age + AgeEffect2Leisure * ((Age*Age)/100) + L;  %Leisure husband

alpha4 = InterceptWelfare + BigCityWelfare * BigCity + SwedishBornWelfare * BornInTheCountry + PrimarySchoolWelfare * PrimarySchool + AgeEffectWelfare18to29 * Age18to29 + AgeEffectWelfare30to34 * Age30to34 + SA; %Welfare

%--------Wife, alfa3-------------

InterceptLeisure = 14.951;
BigCityLeisure = 0.163;
PrimarySchoolLeisure = 2.195;
SecondarySchoolLeisure = 0.552;
AgeEffectLeisure = -0.424;
AgeEffect2Leisure = 0.478;

%FIXED COST
bfcw = 8.507;

%WELFARE - alpha4
InterceptWelfare = 13.826;
BigCityWelfare = -0.094;
SwedishBornWelfare = 3.363;
PrimarySchoolWelfare = -1.702;
AgeEffectWelfare18to29 = -1.857;
AgeEffectWelfare30to34 = -0.647;

%UNOBSERVED HETEROGENEITY
L1 = 14.578; %Leisure
L2 = -7.448;
L3 = -7.130; 

L = L1 * 0.014 + L2 * 0.966 + L3 * 0.02; %Parameter * probability
SA = SA1 * 0.014 + SA2 * 0.966 + SA3 * 0.02; %Parameter * probability
alpha3 = InterceptLeisure + BigCityLeisure * BigCity + PrimarySchoolLeisure * PrimarySchool + SecondarySchoolLeisure * SecondarySchool + AgeEffectLeisure * Age + AgeEffect2Leisure * ((Age*Age)/100) + L;  %Leisure wife
%HÄR SKA ALFA4 VARA

%-------------

bfc = bfcw;

      Psa = 0;
      Dw = 1;
      YearlyWorkIncomeHusband = WageHH * hh;
      C = YearlyWorkIncomeHusband + MonthlyTaxableTransfers * 12 - CalcTaxSwe(YearlyWorkIncomeHusband + MonthlyTaxableTransfers * 12, YearlyWorkIncomeHusband, Age) + MonthlyNonTaxableTransfers * 12; %Netincome husband
      YearlyWorkIncomeWife = WageHW * hw; %TAXABLE TRANSFERS SKA OCKSÅ KÖNSSKILJAS LIKSOM AGE
      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) ... %OBS! Interation fritid och C ska vara gånger 2, se Floods mejl
          + 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 

function Tax = CalcTaxSwe(YearlyIncome, WorkIncome, Age)
global LocalTaxRate
global StateTax1Rate
global StateTax2Rate
%LocalTaxRate = 0.3212; %Average 2017
%StateTax1Rate = 0.25;
%StateTax2Rate = 0.05;
StatTax1Lvl = 430200; 
StatTax2Lvl = 625800;
PBB = 44800; 
IBB = 59300;

      %GRUNDAVDRAG
      if Age < 65
          if floor(YearlyIncome/100)*100<0.99*PBB %Rounded up and down to nearest hundred
              BasicTaxDeductionSwe = ceil(0.423*PBB/100)*100;
          elseif floor(YearlyIncome/100)*100<2.72*PBB
              BasicTaxDeductionSwe = ceil((0.225*PBB+0.2*YearlyIncome)/100)*100;
          elseif floor(YearlyIncome/100)*100<3.11*PBB
              BasicTaxDeductionSwe = ceil(0.77*PBB/100)*100;
          elseif floor(YearlyIncome/100)*100<7.88*PBB
              BasicTaxDeductionSwe = ceil((1.081*PBB-0.1*YearlyIncome)/100)*100;
          else
              BasicTaxDeductionSwe = ceil(0.293*PBB/100)*100;
          end

      else %Older than 64 OTESTAD
          if floor(YearlyIncome/100)*100<0.99*PBB %Rounded up and down to nearest hundred
              BasicTaxDeductionSwe = ceil(0.99*PBB/100)*100;
          elseif floor(YearlyIncome/1)*1<1.11*PBB
              BasicTaxDeductionSwe = ceil(1.11*PBB/100)*100;
          elseif floor(YearlyIncome/10)*10<2.72*PBB
              BasicTaxDeductionSwe = ceil((0.834*PBB+0.249*YearlyIncome)/100)*100;
          elseif floor(YearlyIncome/100)*100<3.77*PBB
              BasicTaxDeductionSwe = ceil(1.511*PBB/100)*100;
          elseif floor(YearlyIncome/100)*100<5.4*PBB
              BasicTaxDeductionSwe = ceil((1.888*PBB-0.1*YearlyIncome)/100)*100;
          elseif floor(YearlyIncome/100)*100<12.43*PBB
              BasicTaxDeductionSwe = ceil((1.834*PBB-0.09*YearlyIncome)/100)*100;
          else
              BasicTaxDeductionSwe = ceil(0.715*PBB/100)*100;
          end
      end

      %EITC
      if Age < 65
          if floor(WorkIncome/100)*100<0.91*PBB %Rounded down and rounded
              EITCSwe = floor((floor(WorkIncome/100)*100-BasicTaxDeductionSwe)*LocalTaxRate);
          elseif floor(WorkIncome/100)*100<2.94*PBB
              EITCSwe = floor((0.91*PBB+0.332*(floor(WorkIncome/100)*100-0.91*PBB)-BasicTaxDeductionSwe)*LocalTaxRate);
          elseif floor(WorkIncome/100)*100<8.08*PBB
              EITCSwe = floor((1.584*PBB+0.111*(floor(WorkIncome/100)*100-2.94*PBB)-BasicTaxDeductionSwe)*LocalTaxRate);
          elseif floor(WorkIncome/100)*100<13.54*PBB
              EITCSwe = floor((2.155*PBB-BasicTaxDeductionSwe)*LocalTaxRate);
          else
              EITCSwe = floor(((2.155*PBB-BasicTaxDeductionSwe)*LocalTaxRate)-0.03*(floor(WorkIncome/100)*100-13.54*PBB)); %Avtappningseffekt
          end

      else %Older than 64
          if floor(WorkIncome/100)*100<100001
              EITCSwe = ceil(0.2*floor(WorkIncome/100)*100);
          elseif floor(WorkIncome/100)*100<300001
              EITCSwe = 15000 + ceil(0.05*floor(WorkIncome/100)*100);
          elseif floor(WorkIncome/100)*100<600001
              EITCSwe = 30000;
          elseif floor(WorkIncome/100)*100<1600001
              EITCSwe = 30000 - ceil(0.03*(floor(WorkIncome/100)*100-600000));
          else
              EITCSwe = 0;
          end
      end

      if EITCSwe < 0 %The EITCSwe cannot be negative!
          EITCSwe = 0;
      end

      %PENSION CONTRIBUTION
      if YearlyIncome > 0.423 * PBB %If the income is below 0.423 PBB no pension contribution!
          if YearlyIncome > 8.07 * IBB %OBS! IBB not PBB
              PensionContribution = round((8.07 * IBB * 0.07)/100)*100; %Only the income above 8.07 IBB is taxed 
          else
              PensionContribution = round((0.07 * YearlyIncome)/100)*100;
          end
      else
          PensionContribution = 0;
      end

      TaxableIncome = floor((YearlyIncome - BasicTaxDeductionSwe)/100) * 100;

      %A = floor(((TaxableIncome - StatTax1Lvl) * StateTax1Rate + (TaxableIncome - StatTax2Lvl) * StateTax2Rate))
      %B = floor(((TaxableIncome - StatTax1Lvl) * StateTax1Rate)/10)*10
      %C = round((TaxableIncome - StatTax1Lvl) * StateTax1Rate + (TaxableIncome - StatTax2Lvl) * StateTax2Rate)

      %CALC TOTAL TAX
      if YearlyIncome > StatTax2Lvl %In order to avoid negative tax rates!
          Tax = floor(TaxableIncome * LocalTaxRate - EITCSwe + floor(((TaxableIncome - StatTax1Lvl) * StateTax1Rate + (TaxableIncome - StatTax2Lvl) * StateTax2Rate)));
      elseif YearlyIncome > StatTax1Lvl
          Tax = floor(TaxableIncome * LocalTaxRate - EITCSwe + floor(((TaxableIncome - StatTax1Lvl) * StateTax1Rate)));
      elseif YearlyIncome > BasicTaxDeductionSwe
          Tax = floor(TaxableIncome * LocalTaxRate - EITCSwe);
      else
          Tax = 0;
      end

      %PensionContribution

      PensionContribution = PensionContribution - Tax;

      if PensionContribution > 0; %The pension contribution is deducted by the tax, but only if there is any tax to deduct it with!
          Tax = Tax + PensionContribution;
      end

      %BasicTaxDeductionSwe
      %EITCSwe
  end

The problem arises when uFamily calls the tax function, since the tax function then think that yearly income = 100 * hh and return an error.

Once again a huge thank you for all your help!

[Walter Roberson]

Really, your function should be using a different husband age than wife age.

Why do you not have something similar to alpha4 for the wife? It is understandable that the male and female constants might be different for welfare effects, but there is a complete expression that is missing for the wife.

If it is the case that only one of the husband or wife can make a particular claim, then it is not clear to an outsider that it must be the husband who makes it -- is there an actual asymmetry built into the tax laws, or is the couple not entitled to use whichever is the most advantageous ?

The husband / wife system with asymmetric coefficients does not account for same-sex marriage, which has been legal in Sweden for about 8 years.

I have been working on the code to use symbolic expressions to allow for the possibility of using an analytic optimizer. Unfortunately parts of it are pretty slow! The individual Swedish tax calculations for husband and wife are not bad at producing symbolic expressions in a reasonable time, but as soon as you add the two together to create your "C", then that is slow, as it works through all the combinations of circumstances. Then with the resulting C being quite long, the calculations of the various uFamily terms is pretty slow.

I have to go out for a while, so it may take some time to finish.

[KGB91]

Yeah, later on I will use micro data and import the variables from real data (like age etc), this is just a test for fixing the optimization problem. Many thanks for your help!!! :)

[Walter Roberson]

Could you confirm that in the section

else %Older than 64 OTESTAD

the lines

          elseif floor(YearlyIncome/1)*1<1.11*PBB
              BasicTaxDeductionSwe = ceil(1.11*PBB/100)*100;
          elseif floor(YearlyIncome/10)*10<2.72*PBB
              BasicTaxDeductionSwe = ceil((0.834*PBB+0.249*YearlyIncome)/100)*100;

Other than those two clauses, everything else in that block and the main "if" use floor(YearlyIncome/100)*100 . The 1.11 line is the only place that you quantize to the nearest kronar, and the 2.72 line is the only place that you quantize to the nearest 10 kronar. It seems especially odd to quantize the income to 1 or 10 kronar but to quantize the deduction to the nearest 100 kronar.

(This is a section that turns out to be important to rewrite for efficiency, so I am reformulating how the calculations are done, and want to be sure that the current code is deliberate rather than a typo.)

[KGB91]

Is should be correct, yes. I will double check it now and return to you asp, but the Swedish tax system is a bit strange sometimes.

[KGB91]

You're right, there is a mistake in the code. IBB should also be 61 500 and not the value in the code (which is last year's value). I'll update the code now.

[KGB91]

Those two things (and the fact that IBB should be 61 500) should be the errors in the code, when I double check with the documents and the tax board's calculations. Thanks! :)

[Walter Roberson]

My symbolic work is still unacceptably slow :(

When I plot the C results, they are nearly planar, except for a small change around 180

[KGB91]

How slow? Before I came here I made a manual code that takes about 10 min per family:

max_i = 0;
max_j = 0;
max_value = 0
for i =1:3500
   for j = 1:3500
      new_value = CalcUFamily(i,j);
      if new_value > max_value;
         max_value = new_value;
         max_i = i;
         max_j = j;
      end
  end
end
max_i
max_j

When finished, I will run the code for about 3000-4000 families and five times (for each tax system). So everything above one minute per family would be a bit slow, but anything under that is fine. I can let the machine run during the night. The upper bound can be lowered to 3000 h or even 2500 h if necessary - even if 3500 h is preferred.

[Walter Roberson]

My symbolic approach is taking more than an hour just to calculate what the function is, and then takes a while to calculate the result of the function for any one pair of values. I have not started algebraic optimization yet.

[KGB91]

That's a while. :(