function [ys,params,check] = EnvironmentA_steadystate(ys,exo,M_,options_)
% 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] 0 if steady state computation worked and to
%                        1 of not (allows to impose restrictions on parameters)
   
%% Step 0: initialize indicator and set options for nuumerical solver
check = 0;
options = optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);

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

% The epsilon parameter is redefined
epssilon = (1-siggma) / siggma;

%% Step 2: Check parameter restrictions
%if ettaC*ettaL<1 % parameter violates restriction (here it is artifical)
%    check=1; %set failure indicator
%    return;  %return without updating steady states
%end
    
%% Step 3: Enter model equations here
% Legend:
    % x(1): N_y, x(2): K_y, x(3): N_ef, x(4): K_ef, x(5): N_er, x(6): K_er,
    % x(7): K_g

Model_wo_Tax = @(x) [(1/rhho - (1-deltta_y)) - (betta * x(1)^alppha * x(2)^(betta-1) * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta));
                 (1/rhho - (1-deltta_ef)) - ((1-thetta) * (xxi * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (xxi + etta * (x(3)^thetta * x(4)^(1-thetta))^epssilon * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon))^((-1/epssilon) - 1)) * x(3)^(thetta) * x(4)^(-thetta));
                 (1/rhho - (1-deltta_er)) - (kapa * (nuu * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (nuu + xxi * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^epssilon * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^((-1/epssilon) - 1)) * x(5)^iotta * x(6)^(kapa-1) * x(7)^(1-iotta-kapa));
                 (1/rhho - (1-deltta_g) - (tauu * (x(3)^thetta * x(4)^(1-thetta))) / x(7)) - ((1-iotta-kapa) * (nuu * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (nuu + xxi * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^epssilon * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^((-1/epssilon) - 1)) * x(5)^iotta * x(6)^kapa * x(7)^(-iotta-kapa));
                 (x(1)^chhi) - (((x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta)) - deltta_y * x(2) - deltta_ef * x(4) - deltta_er * x(6) - deltta_g * x(7))^(-qq) * (alppha * x(1)^(alppha-1) * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta)));
                 (x(3)^ommega) - (((x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta)) - deltta_y * x(2) - deltta_ef * x(4) - deltta_er * x(6) - deltta_g * x(7))^(-qq) * (thetta * (xxi * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (xxi + etta * (x(3)^thetta * x(4)^(1-thetta))^epssilon * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon))^((-1/epssilon) - 1)) * x(3)^(thetta-1) * x(4)^(1-thetta)));
                 (x(5)^pssi) - (((x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta)) - deltta_y * x(2) - deltta_ef * x(4) - deltta_er * x(6) - deltta_g * x(7))^(-qq) * (iotta * (nuu * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (nuu + xxi * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^epssilon * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^((-1/epssilon) - 1)) * x(5)^(iotta-1) * x(6)^kapa * x(7)^(1-iotta-kapa)))];

x0 = [ones(7,1)];

% options = optimoptions('fsolve','Display','iter');

% Computing the steady state
[x, fval] = fsolve(Model_wo_Tax,x0,options) 

% Assign the steady state values to the variables:
N_y = x(1);
K_y = x(2);
N_ef = x(3);
K_ef = x(4);
N_er = x(5);
K_er = x(6);
K_g = x(7);
Y = x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta);
E = (etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon);
EF = x(3)^thetta * x(4)^(1-thetta);
ER = x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa);
W_y = alppha * x(1)^(alppha-1) * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta);
P_e = (1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta);
P_ef = xxi * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (xxi + etta * (x(3)^thetta * x(4)^(1-thetta))^epssilon * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon))^((-1/epssilon) - 1) - tauu;
P_er = nuu * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (nuu + xxi * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^epssilon * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^((-1/epssilon) - 1);
W_ef = thetta * (xxi * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (xxi + etta * (x(3)^thetta * x(4)^(1-thetta))^epssilon * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon))^((-1/epssilon) - 1)) * x(3)^(thetta-1) * x(4)^(1-thetta);
W_er = iotta * (nuu * ((1-alppha-betta) * x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(-alppha-betta)) * (nuu + xxi * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^epssilon * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^((-1/epssilon) - 1)) * x(5)^(iotta-1) * x(6)^kapa * x(7)^(1-iotta-kapa);
C = (x(1)^alppha * x(2)^betta * ((etta * (x(5)^iotta * x(6)^kapa * x(7)^(1-iotta-kapa))^(-epssilon) + xxi * (x(3)^thetta * x(4)^(1-thetta))^(-epssilon))^(-1/epssilon))^(1-alppha-betta)) - deltta_y * x(2) - deltta_ef * x(4) - deltta_er * x(6) - deltta_g * x(7);
Muu = (tauu * (x(3)^thetta * x(4)^(1-thetta))) / x(7);
Rs_y = 1/rhho;
Rs_ef = 1/rhho;
Rs_er = 1/rhho;
Rs_g = 1/rhho;
rs_y = Rs_y - (1-deltta_y);
rs_ef = Rs_ef - (1-deltta_ef);
rs_er = Rs_er - (1-deltta_er);
rs_g = Rs_g - (1-deltta_g);
y_obs = 0;
c_obs = 0;
e_obs = 0;
ef_obs = 0;
er_obs = 0;

%% Step 4: Update parameters and variables
params=NaN(M_.param_nbr,1);
for iter = 1:M_.param_nbr %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

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