% --- Steady state computation --- %
% 
function [ys,params,check] = effort_ramsey_steadystate(ys,exo,M_,options_)
% 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 © 2013-2024 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) ');']);
    eval([ 'params.' paramname '= M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;



%% THIS BLOCK IS MODEL SPECIFIC.
%%
%% Here the user has to define the steady state.
%%

Y          = Y_SS;
n          = n_SS;    
q          = q_SS;  
theta      = theta_SS;   
u          = u_SS;  
v          = theta*u;
s          = s_SS;
mpn        = mpn_SS;
upsilon    = upsilon_SS;

Psiw       = 0;
Psip       = 0; 
Omegaw     = 1;
Omegap     = 1; 
ZI         = 1;
Zr         = 1;
ZR         = 1;
ZL         = 1;
ZB         = 1;

pinf       = 1;              
uk         = 1;  
pk         = 1;
betat      = betta;
etat       = eta;
R          = pinf/betta; 
rk         = (1-betta*(1-deltak))/betta;
aprime     = gamma1;
a          = 0;
ks         = alfa*s*Y/rk;
K          = ks/uk;
I       = deltak*K/ZI;


% options    =   optimset('Display','iter');
% [XX_ss,fval,exitf] = fsolve(@f_effort,[mrsh0 wh0 b0 bigphi0],options,params);
% mrsh    = XX_ss(1);
% wh      = XX_ss(2);
% b       = XX_ss(3);
% bigphi  = XX_ss(4);

% options    =   optimset('Display','iter');
% [XX_ss,fval,exitf] = fsolve(@f_effort_ramsey,[mrsh0 wh0],options,params);
% mrsh    = XX_ss(1);
% wh      = XX_ss(2);

mrsh = 0.887115;
wh = 1.10340496184;


G          = Gy*Y;
C          = Y - I - c*v - G;
biglambda  = 1/((1-lambdac)*C);
h          = ( (sigmae/(1+sigmah+sigmae))*mrsh*biglambda )^(1/(1+sigmah));
eff        = e0*h^(sigmah/(1+sigmae));
gprime     = (1+sigmah+sigmae)/sigmae*h^sigmah;
nh         = n*h;           
prod       = Y/(n*h);
ls         = eff*n*h;                 
A          = Y/(ls^(1-alfa)*ks^alfa);
w          = wh/h;
sL         = w*n/Y;
epsilont   = epsilon;
epsilonor  = (s_SS-1)*epsilon;

Y_dev = 0;
C_dev = 0;
v_dev = 0;
h_dev = 0;
w_dev = 0;
K_dev = 0;
s_dev = 0;
sL_dev = 0;
nh_dev = 0;
prod_dev = 0;
n_dev = 0; 
u_dev = 0;
eff_dev =0;
pinf_dev = 0;
R_dev = 0;
I_dev = 0;
uk_dev = 0;

n_n = n;
theta_n = theta;
q_n = q;
Y_n = Y; 
C_n = C;
h_n = h;
eff_n = eff;
ls_n = ls;
s_n = s;
R_n = pinf/betta;
v_n = v;
u_n = u;
prod_n = prod;
w_n = w;
betat_n = betat;
biglambda_n = biglambda;
I_n = I;
rk_n = rk;
K_n = K;
uk_n = uk;
ks_n = ks;
pk_n = pk;
a_n = 0;
aprime_n = gamma1;
mpn_n = (1-alfa)*Y_n/n_n; 
upsilon_n = upsilon;
mrsh_n = mrsh;
gprime_n = gprime;
ygap  =log(Y)-log(Y_n); 
pi_ann = 4*log(pinf);

dy     = ctrend;
dc     = ctrend;
di     = ctrend ;
dw     = ctrend ;
du     = un;
piobs  = pinf_dev  ;
robs   = conster;
dnh    = consth;



%%
%% END OF THE MODEL SPECIFIC BLOCK.


%% DO NOT CHANGE THIS PART.
%%
%% 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

