% This one tries to compute steady state with only for the variables that do 
% not have analytic closed form

function F = num_ss_Al(x)%, fixed)


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



pssi = 1.52;
gama = 0.05;
theta = 1;
betta = 0.99;
phi_k = 45.5;
phi_p = 0.47;
zzeta = 41.6;
eta_F = 0.088;
eta_I = 0.6;
delta = 0.05/4;
alpha_k = 1/3;
a_I = 0.67;
eps_p = 5;
eps_L = 8;
omega = 0.4;
mu = 0.75;
xi_F = 0.71;
xi_I = 0;
rho_a = 0.92;
g_bar = 0.09;
kappa_F = 0.21;
kappa_I = 0;
rho_int_cost = 0.87;
tau_payroll = 0.25;
pi_ee = 0;
i_ss = (1+pi_ee)/betta - 1;

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_ss;
SDF_t = betta;


F(1) = q_I_t - ((un_t/(eta_I*(L_I_t/q_I_t)))^(omega));

F(2) = zzeta*un_t - (pr_F_t*V_I_h_t/((1/(c_t - pssi*c_t*(L_I_t+L_F_t)^(1+theta)/(1+theta)))*(1-(pssi*gama*(L_I_t+L_F_t)^(1+theta))/((1+theta)*(1-betta*(1-gama))))) ...
	+ pr_I_t*(V_F_h_t/((1/(c_t - pssi*c_t*(L_I_t+L_F_t)^(1+theta)/(1+theta)))*(1-(pssi*gama*(L_I_t+L_F_t)^(1+theta))/((1+theta)*(1-betta*(1-gama)))))));

F(3) = V_F_h_t - (w_F_t - pssi*c_t*(L_I_t+L_F_t)^(theta)*(1/(c_t - pssi*c_t*(L_I_t + L_F_t)^(1+theta)/(1+theta))))/(1-(1-eta_F)*SDF_t);

F(4) = V_I_h_t - ((w_I_t - pssi*c_t*(L_I_t+L_F_t)^(theta)*(1/(c_t - pssi*c_t*(L_I_t + L_F_t)^(1+theta)/(1+theta)))) / (1-(1-eta_I)*SDF_t));

F(5) = L_I_t + L_F_t + un_t + O_t - 1;

F(6) = (a_F_t*L_F_t) - ((p_F_t/(((1-alpha_k)*mc_t*((a_t/ups_p_t)*k_t^(alpha_k)*y_L_t^(1-alpha_k)))/y_L_t))^(-eps_L)*y_L_t);

F(7) = y_L_t - (((a_I*L_I_t)^((eps_L-1)/eps_L)+(a_F_t*L_F_t)^((eps_L-1)/eps_L))^(eps_L/(eps_L-1)));

F(8) = r_t*k_t - (alpha_k*mc_t*((a_t/ups_p_t)*k_t^(alpha_k)*y_L_t^(1-alpha_k)));

F(9) = (a_I*L_I_t) - ((p_I_t/(((1-alpha_k)*mc_t*((a_t/ups_p_t)*k_t^(alpha_k)*y_L_t^(1-alpha_k)))/y_L_t))^(-eps_L)*y_L_t);

F(10) = V_I_f_t - ((p_I_t*a_I - (1+kappa_I*i_l_t)*w_I_t)/(1-(1-eta_I)*SDF_t));

F(11) = V_F_f_t - ((p_F_t*a_F_t - (1+kappa_F*i_l_t)*w_F_t)/(1-(1-eta_F)*SDF_t));

F(12) = q_I_t*V_I_f_t - (xi_I);

F(13) = pr_F_t - ((un_t/eta_F*(L_F_t/q_F_t))^(omega-1));

F(14) = q_F_t - ((un_t/eta_F*(L_F_t/q_F_t))^(omega));

F(15) = (1-mu)*V_F_h_t - (mu*V_F_f_t);

F(16) = q_F_t*V_F_f_t - (xi_F);

F(17) = (1-g_t)*((a_t/ups_p_t)*k_t^(alpha_k)*y_L_t^(1-alpha_k)) - (c_t + (delta*k_t) + xi_F*(eta_F*L_F_t/q_F_t) + xi_I*(eta_I*L_I_t/q_I_t));

F(18) = (1-mu)*V_I_h_t - (mu*V_I_f_t);

F(19) = pr_I_t - ((un_t/(eta_I*(L_I_t/q_I_t)))^(omega-1));



