function [ys,params,check] = baseline_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = NK_baseline_steadystate(ys,exo,M_,options_)
% 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
%   - M_        [sTucture] Dynare model sTucture
%   - options   [sTucture] Dynare options sTucture
%
% 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 resTictions on parameters)

% Copyright (C) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redisTibute 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 disTibuted 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 <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

options=optimset(); % set options for numerical solver

O = 0;
v_m = 0;
r_b = 1/beta_sa-1;
r_a = r_b;
r_k = r_a+delta;
pie = 0;
i = (1+r_b)*(1+pie)-1;
m = (epsilon-1)/epsilon;
lambda_tt1 = beta_sa;
Ty = 0.06;
Bgy  = -4*0.23;
Gy =  r_b*Bgy+tau*((1-alph)*m+(1-omega)*(1-m))-Ty;
Ty_sp = Lambda_T/Lambda*Ty;
Ty_sa = 1/(1-Lambda)*Ty-Lambda/(1-Lambda)*Ty_sp;
by_sa = -Bgy/(1-Lambda);
s = 1/(1-Lambda);
Ky = alph*m/r_k;
w = (1-alph)*m^(1/(1-alph))*(r_k/alph)^(-alph/(1-alph));
Ny = (1-alph)*m/w;
ky_sa = Ky/(1-Lambda);
Piy = (1-m);
qy = omega*Piy/r_a;
ay = ky_sa+s*qy;
dy = -r_a*ay;

fun = @(x)paramfun(x, tau, w, omega, m, Ty_sp, Ty_sa, r_b, by_sa,dy, xi_1, xi_2, Ky, Ny, alph, psi_l, nu, Lambda);
x0 = [0.3 0.3 0.3 0.3 0.3];
x = fsolve(fun,x0);

N  = x(5);
cy_sp = x(1);
ly_sp = x(2);

cy_sa = x(3);
ly_sa = x(4);
Y = (Ky/Ny)^alph*N;
K = Ky*Y;
l_sa = ly_sa*Y;
l_sp = ly_sp*Y;
c_sp = cy_sp*Y;
c_sa = cy_sa*Y;
k_sa = ky_sa*Y;
F_sp = (1-omega)*(1-m)*Y;
F_sa = (1-omega)*(1-m)*Y;
Pi  = (1-m)*Y;
T = Ty*Y;
T_sp =  Ty_sp*Y;
T_sa = Ty_sa*Y;
B_g = Bgy*Y;
G  = Gy*Y;
b_sa = by_sa*Y;
a = ay*Y;
d = -r_a*a;
q = qy*Y;
X = xi_1*abs(d)^(xi_2);
I = delta*K;
C = Y-G-I-(1-Lambda)*X;
eta = 1+sign(d)*xi_1*xi_2*abs(d)^(xi_2-1);
G_bar = G;
B_g_bar = B_g;
T_bar = T;







%% 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
