function [ys, params, check] = alber_urr_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)

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

% Steady state values
x0 = ones(19,1);


options = optimset('Display', 'iter', 'TolX', 1e-3, 'TolFun', 1e-4,...
	'MaxFunEvals',10,'MaxIter',1e7);
	% 'MaxFunEvals',1e7,'MaxIter',1e7);
[x1 Fval] = fsolve(@(x) num_ss_Al(x), x0, options)


w_F_t = x1(1);
V_F_h_t = x1(2);
V_I_h_t = x1(3);
c_t = x1(4);
k_t = x1(5);
O_t = x1(6);
un_t = x1(7);
L_I_t = x1(8);
L_F_t = x1(9);
p_I_t = x1(10);
y_L_t = x1(11);
p_F_t = x1(12);
w_I_t = x1(13);
V_F_f_t = x1(14);
V_I_f_t = x1(15);
pr_F_t = x1(16);
q_F_t = x1(17);
pr_I_t = x1(18);
q_I_t = x1(19);

i_t = i_ss;
pi_t = pi_ee;
r_t = 1/betta+delta-1;
pi_opt_t = (((1+pi_t)^(1-eps_p)-phi_p)/(1-phi_p))^(1/(1-eps_p))-1;
aux_mc = (((1+pi_t)*(1-phi_p*(1+pi_t)^(eps_p-1)*betta))/(1-phi_p*(1+pi_t)^(eps_p)*betta))^(-1);
mc_t = ((eps_p-1)/eps_p)*(1+pi_opt_t)*aux_mc;
ups_p_t = ((1-phi_p)*(1+pi_opt_t)^(-eps_p)*(1+pi_t)^(eps_p))/(1-phi_p*(1+pi_t)^(eps_p));

a_t = 1;
a_F_t = 1;
g_t = g_bar;
i_l_t = i_t;
SDF_t = betta;

L_t = L_I_t + L_F_t;
Y_t = ((a_t/ups_p_t)*k_t^(alpha_k)*y_L_t^(1-alpha_k));

aux_p_1_t = mc_t*Y_t/(1-phi_p*SDF_t*(1+pi_t)^(eps_p));

aux_p_2_t = Y_t/(1-phi_p*SDF_t*(1+pi_t)^(eps_p-1));
Inv_t = delta*k_t;
X_t = c_t;
aux_k = 0;
int_cost_t = 0;

u_c_t = (1/(c_t - pssi*c_t*L_t^(1+theta)/(1+theta)));
y_I_t = a_I*L_I_t;
y_F_t = a_F_t*L_F_t;
lambda_X_t = -pssi*c_t*L_t^(1+theta)*u_c_t/((1+theta)*(1-betta*(1-gama)));

lambda_c_t = (u_c_t*(1-(pssi*gama*(L_I_t+L_F_t)^(1+theta))/((1+theta)*(1-betta*(1-gama)))));

lambda_F_t = V_F_h_t/lambda_c_t;
lambda_I_t = V_I_h_t/lambda_c_t;
p_L_t = (1-alpha_k)*mc_t*Y_t/y_L_t;

vac_F_t = eta_F*L_F_t/q_F_t;
vac_I_t = eta_I*L_I_t/q_I_t;


% End steady state values

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

