function [ys,params,check] = rbc_det4_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = Gali_2010_calib_target_steadystate(ys,exo,M_,options_)
% uses a numerical solver to get the steady state in order to set the required parameters matching the calibration
% 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 (C) 2016-20 Johannes Pfeifer
%
% This 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.
%
% It 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.
%
% For a copy of the GNU General Public License,
% see <http://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;

%% Enter model equations here

efficiency = exo(strmatch('EfficiencyInnovation',M_.exo_names,'exact'))/(1-rho);
Efficiency = effstar*exp(efficiency);
Output_per_unit_of_Capital=((1/betta-1+delta)/alppha)^(1/(1-pssi));
Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/Efficiency)^pssi-alppha)/(1-alppha))^(1/pssi);
Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;

% Compute steady state share of capital.
ShareOfCapital=alppha/(alppha+(1-alppha)*Labour_per_unit_of_Capital^pssi);

% Compute steady state of the endogenous variables.
Labour=1/(1+Consumption_per_unit_of_Labour/((1-alppha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-pssi)));
Consumption=Consumption_per_unit_of_Labour*Labour;
Capital=Labour/Labour_per_unit_of_Capital;
Output=Output_per_unit_of_Capital*Capital;
ExpectedTerm=betta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption)*(alppha*((Output/Capital)^(1-pssi))+1-delta);

%% 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) ') =  ' M_.endo_names{ii} ';']); % be careful that this is written correctly (different if model is already linearized as then all st. st. values are 0)
end
