Fsolve complex results make no sense

13 次查看(过去 30 天)
N/A
N/A 2020-3-4
回答: Alex Sha 2020-3-5
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

回答(2 个)

Walter Roberson
Walter Roberson 2020-3-4
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.

Alex Sha
Alex Sha 2020-3-5
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

类别

Help CenterFile Exchange 中查找有关 Energy Production 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by