function [ys,params,check] = Vnew2b_steadystate(ys,exo,M_,options_)

% function [ys,check] = *model*_steadystate(ys,exo)
% computes the steady state for the associated *model*.mod file using a numerical solver
%
% 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 for the the endogenous variables
%   - check     [scalar] set to 0 if steady state computation worked and to
%                        1 of not (allows to impose restriction on parameters)

check = 0; % set the indicator vaiable (=1 means that the steady state is not found)


mu=0; %dummy value to make parameter name known

%-------------------------------------------------------------------------
% Upload initial values of structural parameters
%-------------------------------------------------------------------------
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end



%-------------------------------------------------------------------------
% Set up the problem (try to use closed-form solutions as much as possible)
%-------------------------------------------------------------------------

Eqs = @(x)...
[...

% 25 unknowns in the system:
% c_e: x(1)         LF_e: x(2)          LS_e: x(3)
% h_e: x(4)         hF_e: x(5)          hS_e: x(6)  
% k_e: x(7)         inv_e: x(8)         b_e: x(9)   

% y: x(10)      rK: x(11)       rH: x(12)       w: x(13)        g: x(14)
% betta: x(15)  kappa: x(16)    zeta: x(17)     d: x(18)    tau_d: x(19)

% c_w: x(20)        LF_w: x(21)          LS_w: x(22)
% h_w: x(23)        hF_w: x(24)          hS_w: x(25)  

% 1. Eq.4 of the model (combined with Eq.3)
%chi = lambda_e * w / (1-LF_e-LS_e)^eta
%chi = lambda_w * w / (1-LF_w-LS_w)^eta
% lambda_e / (1-LF_e-LS_e)^eta = lambda_w / (1-LF_w-LS_w)^eta
% c_e^(-sigm) / (1-LF_e-LS_e)^eta = c_w^(-sigm) / (1-LF_w-LS_w)^eta
x(1)^(-sigm) / (1-x(2)-x(3))^eta - x(20)^(-sigm) / (1-x(21)-x(22))^eta;

% 2. Eq.5 of the model
% (w + d)*LS_e/(rH*hS_e) - (1-theta)/theta
(x(13)+x(18))*x(3) /(x(12)*x(6)) - ((1-theta)/theta);

% 3. Eq.6 of the model
% (w + d)*LS_w/(rH*hS_w) - (1-theta)/theta
(x(13)+x(18))*x(22) /(x(12)*x(25)) - ((1-theta)/theta);

% 4. Eq.7 of the model
% h_e*(1+g) = zeta* (hS_e^theta) * (LS_e^(1-theta))
x(4)*(1 + x(14)) - x(17)*(x(6)^theta) * (x(3)^(1-theta));

% 5. Eq.8 of the model
% h_w*(1+g) = zeta* (hS_w^theta) * (LS_w^(1-theta))
x(23)*(1 + x(14)) - x(17)*(x(25)^theta) * (x(22)^(1-theta));

% 6. Eq.9 of the model
% hF_e + hS_e = h_e
x(5) + x(6) - x(4);

% 7. Eq.10 of the model
% hF_w + hS_w = h_w
x(24) + x(25) - x(23);

% 8. Eq.11 of the model
% betta * theta * h_e/hS_e * (1+g)^(1-sigm) = 1
x(15) * theta * x(4)/x(6) * (1 + x(14))^(1-sigm) - 1;

% 9. Eq.12 of the model
% betta * theta * h_w/hS_w * (1+g)^(1-sigm) = 1
x(15) * theta * x(23)/x(25) * (1 + x(14))^(1-sigm) - 1;

% 10. Eq.13 of the model
% c_e + inv_e + (rST-g)*b_e + d*LS_e = w*LF_e + rH*hF_e + rK*k_e + 0.5*rm
x(1) + x(8) + (rST-x(14))*x(9) + x(18)*x(3) - (x(13)*x(2) + x(12)*x(5) + x(11)*x(7) + 0.5*tau_rm*x(10)); 

% 11. Eq.14 of the model
% c_w + d*LS_w = w*LF_w + rH*hF_w + 0.5*rm
x(20) + x(18)*x(22) - (x(13)*x(21) + x(12)*x(24) + 0.5*tau_rm*x(10)); 

% Alternative for 11. Eq.13-14 combined (economy's market clearing):
% c_e + c_w + inv_e + (rST-g)*b_e + d*(LS_e + LS_w) = y + rm
% x(1) + x(20) + x(8) + (rST-x(14))*x(9) + x(18)*(x(3) + x(22)) - (x(10) + tau_rm*x(10));

% 12. Eq.15 of the model
% k_e = inv_e /(delta+g) 
x(7) - x(8)/(delta+x(14));

% 13. Eq.17 of the model
% kappa = b_e / k_e
x(16) - x(9)/ x(7);

% 14. Eq.20 of the model
%((1+g)^sigm - kappa) / (rK + 1 - delta - kappa*(1+rST)*(1+g)^(-sigm)) = betta
((1+x(14))^sigm - x(16)) / (x(11) + 1 - delta - x(16)*(1+rST)*(1+x(14))^(-sigm)) - x(15);

% 15. Eq.21 of the model
% y = k_e^alph * (hF_e + hF_w)^mu * (LF_e + LF_w)^(1-alph-mu)
x(10) - ((x(7)^alph) * ((x(5) + x(24))^mu) * ((x(2) + x(21))^(1-alph-mu)));

% 16. Eq.22 of the model
% rK = alph * y/k_e
x(11) - alph * x(10)/x(7);

% 17. Eq.23 of the model
% rH = mu * y/(hF_e + hF_w) 
x(12) - mu * x(10)/(x(5) + x(24));

% 18. Eq.24 of the model
% w = (1-alph-mu) * y/(LF_e + LF_w)
x(13) - (1-alph-mu) * x(10)/(x(2) + x(21));

% 19. Eq.29 of the model
% h_e + h_w = 1
x(4) + x(23) - 1;

% 20. Target: stylized fact
% g = g_targ
x(14) - g_targ;

% 21. Target: stylized fact
% b_e = tau_b * y 
x(9) - tau_b * x(10);

% 22. Target: stylized fact
% inv_e = tau_inv * y
x(8) - tau_inv * x(10);

% 23. Targetted value for parameter
% d = tau_d * y/(LS_e + LS_w)
x(18) - x(19) * x(10)/(x(3) + x(22));

% 24. Target: stylized fact
%LF_w + LS_w = 0.5 
x(21) + x(22) - 0.5;
%x(2) + x(3) - 0.5;

% 28. Option 1:
%LF_w = 0.35
x(21) - 0.35;

% 28. Option 2:
%x(19) - tau_d;

...
];



%-------------------------------------------------------------------------
% Run the solver and save the result
%-------------------------------------------------------------------------

% Initial guesses
%x0 = [0.8108 0.3500 0.1500 0.5000 0.2048 0.2952 11.1365 0.2973 0.3806 2.3788 0.0427 1.7422 1.0195 0.0017 0.9865 0.0342 2.2243 1.2663 0.1597 0.8108 0.3500 0.1500 0.5000 0.2048 0.2952];   
x0 = [0.8353 0.2875 0.2125 0.5000 0.2048 0.2952 10.0929 0.2695 0.3449 2.1558 0.0427 1.5790 1.1249 0.0017 0.9865 0.0342 1.9350 0.3373 0.0665 0.8353 0.2875 0.2125 0.5000 0.2048 0.2952];   

%x0 = load('steady_st_init_values.mat');

% Fsolve options
tol      = 1e-5;
Iter     = 100000;
FunEvals = 100000;
options  = optimset('Display', 'on', 'MaxFunEvals', FunEvals, 'MaxIter', Iter, 'TolFun', tol, 'TolX', tol);

% Run fsolve
[ss,fval,exitflag] = fsolve(Eqs, x0, options);
% [ss,fval,exitflag] = fsolve(Eqs, x0.ss, options);

% Exitflag indicates the reason fsolve stopped;  exitflag < 1 means that the solution was not found
if exitflag <1
    check=1;
    return;
end

% Save the solution
save('steady_st_init_values.mat','ss');



%-------------------------------------------------------------------------
% Assign output to model variables
%-------------------------------------------------------------------------

c_e = ss(1)
c_w = ss(20)
LF_e = ss(2)
LF_w = ss(21)
LS_e = ss(3)
LS_w = ss(22)
h_e = ss(4)
h_w = ss(23)
hF_e = ss(5)
hF_w = ss(24)
hS_e = ss(6)
hS_w = ss(25)
k_e = ss(7)
inv_e = ss(8)
b_e = ss(9)
y = ss(10)
rK = ss(11)
rH = ss(12)
w = ss(13)
g = ss(14)
betta = ss(15)
kappa = ss(16)
zeta = ss(17)
d = ss(18)
tau_d = ss(19)
lambda_e = c_e^(-sigm)                                  %Eq.1 of the model
lambda_w = c_w^(-sigm)                                  %Eq.2 of the model
chi = lambda_e * w / (1-LF_e-LS_e)^eta                  %Eq.3 of the model
phi_e = 0                                               %Eq.16 of the model
omega_e = (1 - betta*(1+rST)*(1+g)^(-sigm)) * lambda_e  %Eq.18 of the model
q_e = lambda_e                                          %Eq.19 of the model
rm = tau_rm * y                                         %Target: stylized fact
rm_bar = rm                                             %Eq.25 of the model
AF = 0                                                  %Eq.26 of the model
AS = 0                                                  %Eq.27 of the model
rm_tld = 0                                              %Eq.28 of the model



%-------------------------------------------------------------------------
% Save the structural parameters (the calibrated ones will be updated)
%-------------------------------------------------------------------------
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
