function [ys,params,check] = edsgeFF_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = edsgeFF_steadystate(ys,exo,M_,options_)
% computes the steady state for the edsgeFF
% Inputs:
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output:
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% Copyright © 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
    paramname = M_.param_names{ii};
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%%

global betta R alfag alfab deltab deltag ksi taubar ...
       rhoL theta1 theta2 omega deltax erow rhoy alfayb alfayg Abar d0 d1 d2 
   
tau = exo(strmatch('tau',M_.exo_names,'exact'));
taubar = tau;
  
R = 1/betta;
Rkb=1/betta;
Rkg= Rkb;

options=optimset('MaxFunEvals',10000,'TolX',1e-10,'TolFun',1e-10);

Kb_Lb_g = 30;
Kg_Lg_g = 40;
Lb_g= 0.20;
Lg_g= 0.30;
emis_g = 0.5;

guess = [Kb_Lb_g,Kg_Lg_g,Lb_g,Lg_g,emis_g]; 
ss_sol = fsolve(@findss,guess,options);

Kb_Lb = ss_sol(1);
Kg_Lg = ss_sol(2);
Lb = ss_sol(3);
Lg = ss_sol(4);
emis =ss_sol(5);

Kb = Kb_Lb*Lb;
Kg = Kg_Lg*Lg;
X = (emis+erow)/(1-deltax);
dX = d2*(X)^2+d1*(X)+d0;
Yb = (1-dX) *Abar* Kb^(alfab)*Lb^(1-alfab);
Yg = (1-dX)*Abar*(Kg^(alfag)*Lg^(1-alfag));
Y = ((alfayb^(1/rhoy))*Yb^((rhoy-1)/rhoy)+(alfayg^(1/rhoy))*Yg^((rhoy-1)/rhoy))^(rhoy/(rhoy-1));
pg=(alfayg*Y/Yg)^(1/rhoy);
pb=(alfayb*Y/Yb)^(1/rhoy);
mu = 1-emis/Yb;
Z = Yb*theta1*mu^theta2;
Ig = deltag*Kg;
Ib = deltab*Kb;
Qg = 1;
Qb = 1;
I = Ig+Ib;
C= Y-I-Z;
L = (Lb^(1+rhoL) + Lg^(1+rhoL))^(1/(1+rhoL));
wb = (1-alfab)*(1-dX)*Abar*(Kb^alfab)*(Lb^(-alfab))*(pb-theta1*mu^theta2-taubar*(1-mu));
wg = (1-alfag)*pg*(1-dX)*Abar*(Kg^alfag)*(Lg^(-alfag));

M = betta;
sg = (Qg*Kg)/(Qg*Kg+Qb*Kb);
spreadb = Rkb-R;
spreadg = Rkg-R;
A = Abar;
eA = 0;

    
%% end own model equations

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
    eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
    varname = M_.endo_names{ii};
    eval(['ys(' int2str(ii) ') = ' varname ';']);
end
