Fsolve complex results make no sense

Hey everyone!

I use the function fsolve to solve an equation system for my model. So far, the function works well. But once I use the second version of equation (9) and use instead of just E, fsolve returns a vector of complex numbers, and does not find a solution. In my opinion this does not make sense, as Eand ν are clearly , and Y that is determined through equation (9) is not interlinked with any other equation. Does somebody have an idea where the problem lies?

Best regards

%Call function
fun   = @steadystateEasy;
guess = ones(9,1);
x     = fsolve(fun,guess)

%------------------------------------------------------------------
%------------------------------------------------------------------
function [SS] = steadystateEasy(x)       
%------------------------------------------------------------------
%------------------------------------------------------------------

% Define Variables for fsolve
    P           = x(1);
    p1_hat      = x(2);
    E           = x(3);
    e1          = x(4);
    e2          = x(5);
    e3          = x(6);
    M1          = x(7);
    S           = x(8);
    Y           = x(9);
%------------------------------------------------------------------
% Calibration
%------------------------------------------------------------------
    L           = 1;                % Exogenous Labor
    alpha       = 0.3;              % Total Capital Share       
    nu          = 0.055;            % Total Energy Share       
    beta        = 0.985;            % Discount Factor       
    rho         = -0.058;           % Elasticity of Substitution 
    lambda1     = 0.543;            % Oil Share      
    lambda2     = 0.102;            % Coal Share       
    lambda3     = 0.356;            % Green Resources Share       
    g1          = 1;                % Weight of Oil   
    g2          = 1;                % Weight of Coal   
    g3          = 0;                % Weight of Green Resources   
    tau         = 0.182784;         % (excel) Taxes       
    psi         = 0.0228;           % Emission Forever       
    psi0        = 0.393;            % Emission 1-Period       
    psiL        = 0.2;              % Emission Geometrically      
%------------------------------------------------------------------
% Predetermined Variables
    R           = 338;             % coal reserves
    K           = 39.2;            % capital 
    A           = 55.5;            % technology
    p2_hat      = 0.103351955;     % price coal (excel)
    p3_hat      = 0.6;             % price green (excel)
    d0delta     = psiL + 1;        % calculate the d0 (with s = 0)
    S0          = 596;             % initial carbons tock
    Y0          = 75;              % Initial Output
%------------------------------------------------------------------

% 1) Energy Prices
    SS(1) = P - (    (p1_hat)^(rho/(rho-1))*lambda1^(1/(1-rho)) + ...
                     (p2_hat)^(rho/(rho-1))*lambda2^(1/(1-rho)) + ...
                     (p3_hat)^(rho/(rho-1))*lambda3^(1/(1-rho))) ^ ((rho-1)/rho);

% 2)-3) Energy Demand
    SS(2) =  E-(nu*(A*L^(1-alpha-nu)*K^alpha)/P)^(1/(1-nu));

% 4) Domestic Fuel Demand
    SS(3) = e2 - E*( (P*lambda2)/(p2_hat))^(1/(1-rho));
    SS(4) = e3 - E*( (P*lambda3)/(p3_hat))^(1/(1-rho));

% 5) Oil Demand
    SS(5) = e1 - E*( (P*lambda1)/(p1_hat))^(1/(1-rho));

% 6) Oil Market Clearing
    SS(6) = e1 - (1-beta)*R;
   
% 7)  Law of Motion for the Carbon Stock
    SS(7) = S - ( S0 + d0delta*M1 );

% 8) Emissions
    SS(8) = M1 - (g1*e1+g2*e2+g3*e3);
    
% 9) (PROBLEM!) Gross Output
    SS(9) = Y - A*L^(1-alpha-nu)*K^alpha*E;            % => works fine
  % SS(9) = Y - A*L^(1-alpha-nu)*K^alpha*E^nu;         % => Leads to complex numbers
end

Answer 1

[Walter Roberson]

fsolve does not have any way of putting in bounds constraints, so it has no way to prevent any of the parameters being tested from going negative.

fsolve works numerically, so it must probe locations "past" the zero location on all sides. It could easily end up probing a negative location while looking for a location that makes the value zero.

If you need to constrain E to positive, define E=x(3).^2 and understand that there will then be two potential solutions, and remember to square the 3rd element of the solution vector afterwards.

Answer 2

[Alex Sha]

Hi, Frederic Kluser, with your second case of Eq(9), I get the results below, all seem to be prefect:

p: 1.79143003389998
p1_hat: 1.08736564026132
e: 5.63289811457658
e1: 5.07
e2: 9.65264186961398
e3: 5.96729942681112
m1: 14.722641869614
s: 613.667170243537
y: 183.471688370019

Fevl:
1.17683640610267E-14
9.76996261670138E-15
-2.1316282072803E-14
-6.21724893790088E-15
-1.95399252334028E-14
-4.44089209850063E-15
2.27373675443232E-13
1.95399252334028E-14
1.13686837721616E-13