function [ys,check] = Simulation_15_steadystate(ys,exo)
% function [ys,check] = NK_baseline_steadystate(ys,exo)
% computes the steady state for the NK_baseline.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%
% Output: 
%   - ys        [vector] vector of steady state values fpr the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impos restriction on parameters)

global M_ 
%verbosity = 0;

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


%% Enter model equations here
%We define vectors of parameters for simplicity

sectorn=15;
e=0.01;
SIGMA = 5;
RHO = 0.95;
DELT = 0.97;
LEO=[];
PHI=[];
SIGMAv=[];
BETA=[];
GAMMA=[];
LAMBDA=[];

alphabet =[ "AGR", "MIN","UTI","CON","MAN","WHO","RET","TRA","INF","FIN","BUS","EDU","ART","OTH","GOV"];

for i0 = 1:sectorn
    PHI(end+1) = eval([strcat('PHI_',alphabet(i0))]);
    BETA(end+1) = eval([strcat('BETA_',alphabet(i0))]);
    SIGMAv(end+1) = eval([strcat('SIGMA_',alphabet(i0))]);
 %   KSI_@{sector} = ((1-BETA_@{sector})^((1-BETA_@{sector})*SIGMA_@{sector}/(SIGMA_@{sector}-1))*(BETA_@{sector}^BETA_@{sector}))^(-1);
    eval([strcat('KSI_', alphabet(i0), ' = ((1-BETA_', alphabet(i0), ')^((1-BETA_', alphabet(i0), ')*SIGMA_', alphabet(i0), '/(SIGMA_', alphabet(i0), '-1))*(BETA_', alphabet(i0), '^BETA_', alphabet(i0), '))^(-1);')]);
end


for i0 = 1:sectorn
   for i1 = 1:sectorn
       GAMMA(end+1) =  eval([strcat('GAMMA_',alphabet(i0),alphabet(i1))]);
   end
end
 

LEO=Leontieffbis(BETA,GAMMA);


 for i0 = 1:sectorn
    for i1 = 1:sectorn
       eval([strcat('LEO_', alphabet(i0), alphabet(i1), '= LEO(', int2str(i0), ',', int2str(i1), ');')]);
    end
 end

for i0 = 1:sectorn
    LAMBDA_str = [strcat('LAMBDA_' ,alphabet(i0) ,' =0 ')];
    for i1 = 1:sectorn
        LAMBDA_str = [strcat(LAMBDA_str, '+PHI_' ,alphabet(i1) ,'*LEO_', alphabet(i1), alphabet(i0))];
    end
    eval([strcat(LAMBDA_str, ';')])
end



   for i0 = 1:sectorn
       LAMBDA(end+1) =  eval([strcat('LAMBDA_', alphabet(i0))]);
   end


C = 1;
w=1;
   
     for i0 = 1:sectorn
       eval([strcat('e_' ,alphabet(i0) ,'=1;')]);
       eval([strcat('P_', alphabet(i0), '=1;')]);
       eval([strcat('C_', alphabet(i0), '=PHI_', alphabet(i0), ';')]);
      end
   

   for i0 = 1:sectorn
       eval([strcat('Y_', alphabet(i0),  '= LAMBDA_', alphabet(i0), ';')]);
          0;
   end 
  

  for i0 = 1:sectorn
       eval([strcat('L_', alphabet(i0) ,'=BETA_', alphabet(i0), '*Y_', alphabet(i0), ';')]);
       eval([strcat('M_', alphabet(i0) ,'=(1-BETA_', alphabet(i0), ')*Y_', alphabet(i0), '^(1-1/(SIGMA_' ,alphabet(i0), '));')]);
       eval([strcat('MP_', alphabet(i0) ,'=1-BETA_' ,alphabet(i0), ';')]);
       eval([strcat('Part_' ,alphabet(i0), '=PHI_', alphabet(i0), ';')]);
  end
  
   
   for i0 = 1:sectorn
    for i1 = 1:sectorn
       eval([strcat('M_', alphabet(i0), alphabet(i1), '=GAMMA_' ,alphabet(i0), alphabet(i1), '*Y_' ,alphabet(i1), ';')]);
    end
   end
    
   
    for i0 = 1:sectorn
      eval([strcat('lY_', alphabet(i0), '=log(Y_', alphabet(i0), ');')]);
    end
    
    for i0 = 1:sectorn
       eval([strcat('Phat_', alphabet(i0), '=0;')]);
    end
    
   for i0 = 1:sectorn
    for i1 = 1:sectorn
       eval([strcat('theta_', alphabet(i0), alphabet(i1), '=1;')]);
    end
   end




%% end own model equations

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

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