optimization - using fmin

have a non linear maximization problem in Matlab (Max utility):

A=(((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*
    (1+r)+Epsilon_A_S*S_A-(Gamma/2)*(Epsilon_A_S^2)*(Sigma_i^2)

B=(((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+
   (Epsilon_A_S+Epsilon_A_c)*S_A-(Gamma/2)*(Epsilon_A_S+Epsilon_A_c)^2*(Sigma_i^2)-
   Epsilon_A_c*K

z=(K-S_A)/Sigma_i

Utility_A=-exp(-Gamma*A)*normcdf(z+Gamma*Epsilon_A_S*Sigma_i)
  +-exp(-Gamma*B)*(1-normcdf(z+Gamma*(Epsilon_A_S+Epsilon_A_c)*Sigma_i))

I want to maximize utility function by finding: Epsilon_A_S, Epsilon_A_c

The constraints are:

Epsilon_A_S<=2

Epsilon_A_c>=-1

All the rest are parameters.

How do I do I set this problem in Matlab?

[Rik]

Sent by email:

"Dear Rik, Thank you very much for helping me with fminseach problem, as I'm doing my first steps in matlab. I still having problem for finding this equilibrium (I cannot find a feasible solution while trying to get the results in this paper in figure 1: https://www.scirp.org/journal/PaperInformation.aspx?PaperID=74734 can you please help me with that? kind regards, Yossi"

The advice offered on this page is applicable here as well. I won't invest the time to read an academic paper outside my field just to help you write code. If you write a good description of what you want to implement, you can post your question on this forum. Make sure to include what you tried and what specifically is different from your intended result. Have a read here and here. It will greatly improve your chances of getting an answer.

Answer 1

[Rik]

If you can set a finite lower bound of Epsilon_A_S and a finite upper bound of Epsilon_A_c, you can use fminbnd (you will need to re-write your function a bit so you have one vector as input with one set of bounds). If you don't, you'll need the Optimization Toolbox which has the fmincon function.

[Yossi]

But how could I re-write this to one vector if it non-linear?

[Rik]

I was mistaken: with the fminbnd function you can't, but with the fminsearch function, you can. I don't have the values for your inputs V I S_Initial Call_Price r S_A Gamma Sigma_i and K, so I can't really test the code.

A=@(Epsilon_A_S,Epsilon_A_c) ...
    (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
    Epsilon_A_S*S_A-(Gamma/2)*(Epsilon_A_S^2)*(Sigma_i^2);
B=@(Epsilon_A_S,Epsilon_A_c) ...
    (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
    (Epsilon_A_S+Epsilon_A_c)*S_A-...
    (Gamma/2)*(Epsilon_A_S+Epsilon_A_c)^2*(Sigma_i^2)-...
    Epsilon_A_c*K;
z=(K-S_A)/Sigma_i;
Utility_A=@(Epsilon_A_S,Epsilon_A_c)...
    -exp(-Gamma*A(Epsilon_A_S,Epsilon_A_c))*...
    normcdf(z+Gamma*Epsilon_A_S*Sigma_i)+...
    -exp(-Gamma*B(Epsilon_A_S,Epsilon_A_c))*...
    (1-normcdf(z+Gamma*(Epsilon_A_S+Epsilon_A_c)*Sigma_i));

%Take the negative of the utility. If we minimize that, we maximize
%utitliy. This function is unconstrained, but when either epsilon goes over
%the bound value, the cost function value goes to inf.
minimize_function=@(Epsilon_A_S,Epsilon_A_c)...
    -Utility_A(Epsilon_A_S,Epsilon_A_c)+...
    1/(Epsilon_A_S<=2)-1+...
    1/(Epsilon_A_c>=-1)-1;
intial_Epsilon_A_S=0;intial_Epsilon_A_c=0;
fitted_val=...
    fminsearch(@(val)minimize_function(val(1),val(2)),...
    [intial_Epsilon_A_S,intial_Epsilon_A_c]);
[fitted_Epsilon_A_S,fitted_Epsilon_A_c]=deal(fitted_val(1),fitted_val(2));

[Yossi]

I run this script, wasn't good somehow:

S=10;
H_mu=0;
S_A=S+H_mu;
S_B=S-H_mu;
Sigma_i=1.5;
K=10;
T=1;
r=0.02;
Gamma=0.1;
V=2;I=2;

S_Initial=(S-Gamma*(Sigma_i^2)*(V/I))/(1+r)

z=(K-S)/Sigma_i

Call_Price=((S-Gamma*(Sigma_i^2)*(V/I)-K)*(1-normcdf(z+Gamma*Sigma_i*(V/I)))+normpdf(z+Gamma*Sigma_i*(V/I))*Sigma_i)/(1+r)

% Intial Vals:

Epsilon_A_S=(V/I);

Epsilon_A_c=0;

A=@(Epsilon_A_S,Epsilon_A_c) ...
    (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
    Epsilon_A_S*S_A-(Gamma/2)*(Epsilon_A_S^2)*(Sigma_i^2);
B=@(Epsilon_A_S,Epsilon_A_c) ...
    (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
    (Epsilon_A_S+Epsilon_A_c)*S_A-...
    (Gamma/2)*(Epsilon_A_S+Epsilon_A_c)^2*(Sigma_i^2)-...
    Epsilon_A_c*K;
z=(K-S_A)/Sigma_i;
Utility_A=@(Epsilon_A_S,Epsilon_A_c)...
    -exp(-Gamma*A(Epsilon_A_S,Epsilon_A_c))*...
    normcdf(z+Gamma*Epsilon_A_S*Sigma_i)+...
    -exp(-Gamma*B(Epsilon_A_S,Epsilon_A_c))*...
    (1-normcdf(z+Gamma*(Epsilon_A_S+Epsilon_A_c)*Sigma_i));

%Take the negative of the utility. If we minimize that, we maximize
%utitliy. This function is unconstrained, but when either epsilon goes over
%the bound value, the cost function value goes to inf.
minimize_function=@(Epsilon_A_S,Epsilon_A_c)...
    -Utility_A(Epsilon_A_S,Epsilon_A_c)+...
    1/(Epsilon_A_S<=2)-1+...
    1/(Epsilon_A_c>=-1)-1;
intial_Epsilon_A_S=0;intial_Epsilon_A_c=0;
fitted_val=...
    fminsearch(@(val)minimize_function(val(1),val(2)),...
    [intial_Epsilon_A_S,intial_Epsilon_A_c]);
[fitted_Epsilon_A_S,fitted_Epsilon_A_c]=deal(fitted_val(1),fitted_val(2));

any suggestion why?

[Yossi]

In addition, how can I put the constraints: Epsilon_A_S<=2 Epsilon_A_c>=-1

[Rik]

Have you actually looked at the code I wrote? I even added comments. The constraints are added in the minimize_function and you forgot to edit the initial values. And what do you mean with "not good somehow"?

BTW, next time, select your code and hit the {}Code button. It will make your code much more readable and it will make it easier to copy and paste to Matlab.

[Yossi]

I'm sorry, I had my mistake... why have you made:

1/(Epsilon_A_S<=2)-1
and not?

Epsilon_A_S<=2

[Rik]

Because when Epsilon_A_S<=2 becomes false (so when the solver tries a non-valid solution), 1/(Epsilon_A_S<=2)-1 becomes inf. Because the solver tries to minimize the function value, this means that it will only ever find a non-valid solution when the cost function is inf everywhere. For solutions that do comply with the conditions, this extra part becomes 0, so it has no effect on the rest of the cost function.

Answer 2

[Yossi]

How can I another constraint that

Epsilon_A_S>=-1

[Rik]

Is the pattern clear? You need to add the 1/(test)-1 with your additional constraint. I've added it in the lines below.

minimize_function=@(Epsilon_A_S,Epsilon_A_c)...
    -Utility_A(Epsilon_A_S,Epsilon_A_c)+...
    1/(Epsilon_A_S<=2)-1+...
    1/(Epsilon_A_S>=-1)-1+...
    1/(Epsilon_A_c>=-1)-1;

[Yossi]

I got the pattern, thank you! but now I have another strange thing: it is I added another criteria in a "while" loop but it seems not to be fulfilled (the sum of the quantities). Do you have any idea why the loop stop but the criteria doesn't work? Thanks!

clear all; close all
Cond3=0;error_model=0.0000001;

      S=10;
      H_mu=0.6;
      S_A=S+H_mu;
      S_B=S-H_mu;
      Sigma_i=1.5;
      K=10;
      T=1;
      r=0.02;
      Gamma=0.1;
      V=2;I=2;
      fitted_Epsilon_B_S=(V/I);fitted_Epsilon_B_c=0;
      fitted_Epsilon_A_S=(V/I);fitted_Epsilon_A_c=0;
      %%Picking Reasonable initial candidates
      %Eq. (4)
      S_Initial=(S-Gamma*(Sigma_i^2)*(V/I))/(1+r)
      %Eq. (5)
      z=(K-S)/Sigma_i
      Call_Price=((S-Gamma*(Sigma_i^2)*(V/I)-K)*(1-normcdf(z+Gamma*Sigma_i*(V/I)))+normpdf(z+Gamma*Sigma_i*(V/I))*Sigma_i)/(1+r)
  while  Cond3<2 
      %%Player A
      A=@(Epsilon_A_S,Epsilon_A_c) ...
          (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
          Epsilon_A_S*S_A-(Gamma/2)*(Epsilon_A_S^2)*(Sigma_i^2);
      B=@(Epsilon_A_S,Epsilon_A_c) ...
          (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
          (Epsilon_A_S+Epsilon_A_c)*S_A-...
          (Gamma/2)*(Epsilon_A_S+Epsilon_A_c)^2*(Sigma_i^2)-...
          Epsilon_A_c*K;
      z=(K-S_A)/Sigma_i;
      Utility_A=@(Epsilon_A_S,Epsilon_A_c)...
          -exp(-Gamma*A(Epsilon_A_S,Epsilon_A_c))*...
          normcdf(z+Gamma*Epsilon_A_S*Sigma_i)+...
          -exp(-Gamma*B(Epsilon_A_S,Epsilon_A_c))*...
          (1-normcdf(z+Gamma*(Epsilon_A_S+Epsilon_A_c)*Sigma_i));

      %Take the negative of the utility. If we minimize that, we maximize
      %utitliy. This function is unconstrained, but when either epsilon goes over
      %the bound value, the cost function value goes to inf.
      minimize_function=@(Epsilon_A_S,Epsilon_A_c)...
          -Utility_A(Epsilon_A_S,Epsilon_A_c)+...
          1/(Epsilon_A_S<=2)-1+...
          1/(Epsilon_A_S>=-1)-1+...
          1/(Epsilon_A_c<=1)-1+...
          1/(Epsilon_A_c>=-1)-1;
      intial_Epsilon_A_S=fitted_Epsilon_A_S;intial_Epsilon_A_c=fitted_Epsilon_A_c;
      fitted_val=...
          fminsearch(@(val)minimize_function(val(1),val(2)),...
          [intial_Epsilon_A_S,intial_Epsilon_A_c]);
      [fitted_Epsilon_A_S,fitted_Epsilon_A_c]=deal(fitted_val(1),fitted_val(2));

      %%Player B
      A=@(Epsilon_B_S,Epsilon_B_c) ...
          (((V/I)-Epsilon_B_S)*S_Initial-Epsilon_B_c*Call_Price)*(1+r)+...
          Epsilon_B_S*S_B-(Gamma/2)*(Epsilon_B_S^2)*(Sigma_i^2);
      B=@(Epsilon_B_S,Epsilon_B_c) ...
          (((V/I)-Epsilon_B_S)*S_Initial-Epsilon_B_c*Call_Price)*(1+r)+...
          (Epsilon_B_S+Epsilon_B_c)*S_B-...
          (Gamma/2)*(Epsilon_B_S+Epsilon_B_c)^2*(Sigma_i^2)-...
          Epsilon_B_c*K;
      z=(K-S_B)/Sigma_i;
      Utility_B=@(Epsilon_B_S,Epsilon_B_c)...
          -exp(-Gamma*A(Epsilon_B_S,Epsilon_B_c))*...
          normcdf(z+Gamma*Epsilon_B_S*Sigma_i)+...
          -exp(-Gamma*B(Epsilon_B_S,Epsilon_B_c))*...
          (1-normcdf(z+Gamma*(Epsilon_B_S+Epsilon_B_c)*Sigma_i));

      %Take the negative of the utility. If we minimize that, we maximize
      %utitliy. This function is unconstrained, but when either epsilon goes over
      %the bound value, the cost function value goes to inf.
      minimize_function=@(Epsilon_B_S,Epsilon_B_c)...
          -Utility_B(Epsilon_B_S,Epsilon_B_c)+...
          1/(Epsilon_B_S<=2)-1+...
          1/(Epsilon_B_S>=-1)-1+...
          1/(Epsilon_B_c<=1)-1+...
          1/(Epsilon_B_c>=-1)-1;
      intial_Epsilon_B_S=fitted_Epsilon_B_S;intial_Epsilon_B_c=fitted_Epsilon_B_c;
      fitted_val=...
          fminsearch(@(val)minimize_function(val(1),val(2)),...
          [intial_Epsilon_B_S,intial_Epsilon_B_c]);
      [fitted_Epsilon_B_S,fitted_Epsilon_B_c]=deal(fitted_val(1),fitted_val(2));

      if fitted_Epsilon_A_S+fitted_Epsilon_B_S>V+error_model
          S_Initial=S_Initial+0.0001;
          Cond1=0;
      else
          if fitted_Epsilon_A_S+fitted_Epsilon_B_S<V-error_model
              S_Initial=S_Initial-0.0001;
              Cond1=0;
          end

          Cond1=1;
      end

      if fitted_Epsilon_A_c+fitted_Epsilon_B_c>error_model
          Call_Price=Call_Price+0.0001;
          Cond2=0;        
      else
          if fitted_Epsilon_A_c+fitted_Epsilon_B_c<error_model
              Call_Price=Call_Price-0.0001;
              Cond2=0;

          end
          Cond2=1;

      end
      Cond3=Cond1+Cond2

end

Answer 3

[Yossi]

I've tried another code. I know the answer should be: Call_price=0.6, S_Initial=~9.8. but somehow it does not work

clear all; close all
Cond3=0;error_model=0.0000001;

      S=10;
      H_mu=0.6;
      S_A=S+H_mu;
      S_B=S-H_mu;
      Sigma_i=1.5;
      K=10;
      T=1;
      r=0.02;
      Gamma=0.1;
      V=2;I=2;
      fitted_Epsilon_B_S=(V/I);fitted_Epsilon_B_c=0;
      fitted_Epsilon_A_S=(V/I);fitted_Epsilon_A_c=0;
      %%Picking Reasonable initial candidates
      %Eq. (4)
      S_Initial=(S-Gamma*(Sigma_i^2)*(V/I))/(1+r)
      %Eq. (5)
      z=(K-S)/Sigma_i
      Call_Price=((S-Gamma*(Sigma_i^2)*(V/I)-K)*(1-normcdf(z+Gamma*Sigma_i*(V/I)))+normpdf(z+Gamma*Sigma_i*(V/I))*Sigma_i)/(1+r)
  while  Cond3<2 
      %%Player A
      A=@(Epsilon_A_S,Epsilon_A_c) ...
          (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
          Epsilon_A_S*S_A-(Gamma/2)*(Epsilon_A_S^2)*(Sigma_i^2);
      B=@(Epsilon_A_S,Epsilon_A_c) ...
          (((V/I)-Epsilon_A_S)*S_Initial-Epsilon_A_c*Call_Price)*(1+r)+...
          (Epsilon_A_S+Epsilon_A_c)*S_A-...
          (Gamma/2)*(Epsilon_A_S+Epsilon_A_c)^2*(Sigma_i^2)-...
          Epsilon_A_c*K;
      z=(K-S_A)/Sigma_i;
      Utility_A=@(Epsilon_A_S,Epsilon_A_c)...
          -exp(-Gamma*A(Epsilon_A_S,Epsilon_A_c))*...
          normcdf(z+Gamma*Epsilon_A_S*Sigma_i)+...
          -exp(-Gamma*B(Epsilon_A_S,Epsilon_A_c))*...
          (1-normcdf(z+Gamma*(Epsilon_A_S+Epsilon_A_c)*Sigma_i));

      %Take the negative of the utility. If we minimize that, we maximize
      %utitliy. This function is unconstrained, but when either epsilon goes over
      %the bound value, the cost function value goes to inf.
      minimize_function=@(Epsilon_A_S,Epsilon_A_c)...
          -Utility_A(Epsilon_A_S,Epsilon_A_c)+...
          1/(Epsilon_A_S<=2)-1+...
          1/(Epsilon_A_S>0)-1+...
          1/(Epsilon_A_c<=1)-1+...
          1/(Epsilon_A_c>=-1)-1;
      intial_Epsilon_A_S=fitted_Epsilon_A_S;intial_Epsilon_A_c=fitted_Epsilon_A_c;
      fitted_val=...
          fminsearch(@(val)minimize_function(val(1),val(2)),...
          [intial_Epsilon_A_S,intial_Epsilon_A_c]);
      [fitted_Epsilon_A_S,fitted_Epsilon_A_c]=deal(fitted_val(1),fitted_val(2));

      %%Player B
      A=@(Epsilon_B_S,Epsilon_B_c) ...
          (((V/I)-Epsilon_B_S)*S_Initial-Epsilon_B_c*Call_Price)*(1+r)+...
          Epsilon_B_S*S_B-(Gamma/2)*(Epsilon_B_S^2)*(Sigma_i^2);
      B=@(Epsilon_B_S,Epsilon_B_c) ...
          (((V/I)-Epsilon_B_S)*S_Initial-Epsilon_B_c*Call_Price)*(1+r)+...
          (Epsilon_B_S+Epsilon_B_c)*S_B-...
          (Gamma/2)*(Epsilon_B_S+Epsilon_B_c)^2*(Sigma_i^2)-...
          Epsilon_B_c*K;
      z=(K-S_B)/Sigma_i;
      Utility_B=@(Epsilon_B_S,Epsilon_B_c)...
          -exp(-Gamma*A(Epsilon_B_S,Epsilon_B_c))*...
          normcdf(z+Gamma*Epsilon_B_S*Sigma_i)+...
          -exp(-Gamma*B(Epsilon_B_S,Epsilon_B_c))*...
          (1-normcdf(z+Gamma*(Epsilon_B_S+Epsilon_B_c)*Sigma_i));

      %Take the negative of the utility. If we minimize that, we maximize
      %utitliy. This function is unconstrained, but when either epsilon goes over
      %the bound value, the cost function value goes to inf.
      minimize_function=@(Epsilon_B_S,Epsilon_B_c)...
          -Utility_B(Epsilon_B_S,Epsilon_B_c)+...
          1/(Epsilon_B_S<=2)-1+...
          1/(Epsilon_B_S>0)-1+...
          1/(Epsilon_B_c<=1)-1+...
          1/(Epsilon_B_c>=-1)-1;
      intial_Epsilon_B_S=fitted_Epsilon_B_S;intial_Epsilon_B_c=fitted_Epsilon_B_c;
      fitted_val=...
          fminsearch(@(val)minimize_function(val(1),val(2)),...
          [intial_Epsilon_B_S,intial_Epsilon_B_c]);
      [fitted_Epsilon_B_S,fitted_Epsilon_B_c]=deal(fitted_val(1),fitted_val(2));

      if fitted_Epsilon_A_S+fitted_Epsilon_B_S>V+error_model
          S_Initial=S_Initial+0.0001;
          Cond1=0;
      else
          if fitted_Epsilon_A_S+fitted_Epsilon_B_S<V-error_model
              S_Initial=S_Initial-0.0001;
              Cond1=0;
          end

          Cond1=1;
      end

      if fitted_Epsilon_A_c+fitted_Epsilon_B_c>error_model
          Call_Price=Call_Price+0.0001;
          Cond2=0;        
      else
          if fitted_Epsilon_A_c+fitted_Epsilon_B_c<error_model
              Call_Price=Call_Price-0.0001;
              Cond2=0;

          end
          Cond2=1;

      end
      Cond3=Cond1+Cond2
      S_Initial
      Call_Price
  end

[Rik]

You shouldn't use the answer field for comments. And about your code, I think it would help to add an else statement, change this

if fitted_Epsilon_A_S+fitted_Epsilon_B_S<V-error_model
    S_Initial=S_Initial-0.0001;
    Cond1=0;
end
Cond1=1;

to this

if fitted_Epsilon_A_S+fitted_Epsilon_B_S<V-error_model
    S_Initial=S_Initial-0.0001;
    Cond1=0;
else
    Cond1=1;
end

You should really take note of the m-lint warnings, they are very useful in picking up this type of bug. If this doesn't solve your problem, use the debugger to go through your code step by step to see where unexpected things start happening.

[Yossi]

I have a question to your original code: How can I see the value of

Utility_A

[Rik]

It is an anonymous function, so you will have to input the results of the fit to get a value:

Utility_A (fitted_val(1),fitted_val(2))