%Calculates 
function  [ys,params,check] = EDSGE_GK11_steadystate(ys,exo,M_,options_)
%global M_ 
%% DO NOT CHANGE THIS PART.
%%
%% Here we load the values of the deep parameters in a loop.
%%
NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  %varname = deblank(M_.endo_names(ii,:));
  varname = M_.endo_names{ii};
  eval([varname '=ys(' int2str(ii) ');']);
end

NumberOfParameters = M_.param_nbr;                            % Number of deep parameters.
for ii = 1:NumberOfParameters                                  % Loop...
  %paramname = deblank(M_.param_names(ii,:));                   %    Get the name of parameter i. 
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);         %    Get the value of parameter i.
  eval([ 'params.' paramname '= M_.params(' int2str(ii) ');']);
end                                                           % End of the loop.  
check = 0;
%%
%%
%global betta gama R ksi zetta kappa d0 d1 d2 alfag alfab deltab deltag taub taug taubar ...
%       rhoL theta1 theta2 deltax omega erow rhoy alfayb alfayg Abar

taubar = tau;




%Creating new structure for submitted parameters for findssFF_TAX
params_f    =   params;


error_flag=0;
options=optimset('MaxFunEvals',10000,'TolX',1e-20,'TolFun',1e-20);
guess = [Kb_g,Kg_g,Lb_g,E_g];

try
[ss_sol, fval, exitflag]= fsolve(@findssFF_TAX,guess,options,params_f,taubar);
if exitflag<1
   error_flag=1;
end
catch
    error_flag=1;
    ss_sol=guess;
end

Kb = ss_sol(1);
Kg = ss_sol(2);
Lb = ss_sol(3);
E = ss_sol(4);



R   = 1/betta;
Lambda = betta;
X   = (E+erow)/(1-deltax);
dX  = d2*(X)^2+d1*(X)+d0;
Lg = (1 / 3^(1 + rhoL) - Lb^(1 + rhoL))^(1 / (1 + rhoL));
Yb  = (1-dX)*Abar*Kb^(alfab)*Lb^(1-alfab);
Yg  = (1-dX)*Abar*Kg^(alfag)*Lg^(1-alfag);
Y   = ((1-alfayg)^(1/rhoy)*Yb^((rhoy-1)/rhoy)+(alfayg)^(1/rhoy)*Yg^((rhoy-1)/rhoy))^(rhoy/(rhoy-1));
pb  =((1-alfayg)*Y/Yb)^(1/rhoy);
pg  =(alfayg*Y/Yg)^(1/rhoy);
Ig  = deltag*Kg;
Ib  = deltab*Kb;
I   = Ib+Ig;
C   = Y-I;
Qg  = 1;
Qb  = 1;
Rkg = (pg*(1-dX)*alfag*Abar*(Kg/Lg)^(alfag-1) + Qg*(1-deltag))/Qg;
Rkb = (1-dX)*Abar*(Kb/Lb)^(alfab-1)*alfab*(pb - taubar*gamaE)/Qb + (1-deltab);
w  = (1-alfab)*Abar*(1-dX)*(Kb^alfab)*(Lb^(-alfab))*(pb-taubar*gamaE);

p   = (Rkb-(1+taub)*R);
varalfa = (betta-gama)/(gama*betta*p+betta*zetta);
varphi  = ((1-gama)*(betta*p*varalfa+1))/(1-gama-betta*gama*p*varalfa);
nu  = betta*(1-gama+gama*varphi)*R;
chib    = betta*(1-gama+gama*varphi)*(Rkb-(1+taub)*R);
chig    = betta*(1-gama+gama*varphi)*(Rkg-(1+taug)*R);
N   = (gama*((Rkg-(1+taug)*R)*Qg*Kg+(Rkb-(1+taub)*R)*Qb*Kb)+zetta*(Qg*Kg+Qb*Kb))/(1-gama*R);
D   = (1+taub)*Qb*Kb+(1+taug)*Qg*Kg-N;

LbAdj = 0;
LgAdj = 0;

A   = Abar;
eA  = 0;


%%
%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
%% Here we define the steady state values of the endogenous variables of
%% the model.
%%
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;                % Number of endogenous variables.
%ys = zeros(NumberOfEndogenousVariables,1);                     % Initialization of ys (steady state).
for ii = 1:NumberOfEndogenousVariables                         % Loop...
  %varname = deblank(M_.endo_names(ii,:));                     %    Get the name of endogenous variable i.                     
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);                %    Get the steady state value of this variable.
end

