function [ys,params,check] = full_calibrate_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_        [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) 2013-2020 Dynare Team
%
% This file is part of Dynare.
%
% Dynare 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.
%
% Dynare 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.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://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

sol=fsolve(@(x)steadystate(x,ssigma,ppsi,oomegac,tthetac,ggamma,pphi,oomegal,tthetal,xxil,AH,aalphah,AN,aalphan),zeros(1,13));
C = sol(1);
CT = sol(2);
CN = sol(3);
CH = sol(4);
CF = sol(5);
ELL = sol(6);
ELLH = sol(7);
ELLN = sol(8);
YH = log(AH) + aalphah*ELLH;
tildeCT = oomegac^(1/tthetac)*exp(CT)^((tthetac-1)/tthetac)/exp(C)^((tthetac-1)/tthetac);
tildeCN = (1-oomegac)^(1/tthetac)*exp(CN)^((tthetac-1)/tthetac)/exp(C)^((tthetac-1)/tthetac);
tildeCH = ggamma^(1/pphi)*exp(CH)^((pphi-1)/pphi)/exp(CT)^((pphi-1)/pphi);
tildeCF = (1-ggamma)^(1/pphi)*exp(CF)^((pphi-1)/pphi)/exp(CT)^((pphi-1)/pphi);
tildeLH = oomegal^(1/tthetal)*exp(ELLH)^((tthetal-1)/tthetal)/exp(ELL)^((tthetal-1)/tthetal);
tildeLN = (1-oomegal)^(1/tthetal)*exp(ELLN)^((tthetal-1)/tthetal)/exp(ELL)^((tthetal-1)/tthetal);
UF = exp(C)^(-ssigma+1/tthetac)*exp(CT)^(-1/tthetac+1/pphi)*exp(CF)^(-1/pphi)*oomegac^(1/tthetac)*(1-ggamma)^(1/pphi);
c=0;
ct=0;
cn=0;
ch=0;
cf=0;
ell=0;
ellh=0;
elln=0;
nfa=0;
bs=Bs;
is=nnu;
rho=0;
pfs=0;
phs=0;
yh=0;
yn=0;
pin=0;
piwh=0;
piwn=0;
ppn=0;
pnpt=0;
pte=0;
wnpn=0;
whe=0;
i=nnu;
rer=0;


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